Section 2 Scenario Analysis: WCS Priority Landscapes

# turn off the s2 processing 
## https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data
sf::sf_use_s2(FALSE)

See the USFS Wildfire Crisis Strategy (WCS) document for background information.

2.1 Scenarios Considered

There are three scenarios considered in this analysis of the constraints to mechanical forest health and fuel reduction treatments:

# mechanical mgmt allowed if:
n_sc_temp <- 3
data.frame(
  scenario = 1:3
  , cover = rep("Forest or Shrubland",n_sc_temp)
  , protected = rep("Not Protected",n_sc_temp)
  , slope = c("<40%","<60%","<60%")
  , road = c("<1,000ft","<2,000ft","<2,000ft")
  , riparian = c(">100ft",">100ft",">50ft")
  , admin = c("No Designation","No Designation","Any Designation")
) %>% 
  tidyr::pivot_longer(
    cols = -c(scenario)
  ) %>% 
  tidyr::pivot_wider(
    names_from = scenario
    , values_from = value
    , names_prefix = "scenario_"
  ) %>% 
  dplyr::mutate(
    name = factor(
      name
      , levels = c(
        "cover"
        , "protected"
        , "slope"
        , "road"
        , "riparian"
        , "admin"
      )
      , labels = c(
        "NLCD Cover Type"
        , "Protected or<br>IRA Status"
        , "Slope"
        , "Distance to<br>Nearest Road"
        , "Riparian<br>Buffer"
        , "Administrative<br>Designation"
      )
      , ordered = T
    )
  ) %>% 
  dplyr::arrange(name) %>% 
  dplyr::mutate(
    lvl = paste0("L",dplyr::row_number()-1,":")
  ) %>% 
  dplyr::relocate(lvl) %>% 
# make table
kableExtra::kable(
    caption = "Scenarios Considered for Determining Mechanical Treatment Operability<br>Mechanical treatment allowed if:"
    , col.names = c(
      ""
      , ""
      , "Scenario 1<br>Most Constrained"
      , "Scenario 2<br>Moderate"
      , "Scenario 3<br>Least Constrained"
    )
    , escape = F
  ) %>%  
  kable_classic(full_width=T) %>% 
  kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE)
Table 2.1: Scenarios Considered for Determining Mechanical Treatment Operability
Mechanical treatment allowed if:
Scenario 1
Most Constrained
Scenario 2
Moderate
Scenario 3
Least Constrained
L0: NLCD Cover Type Forest or Shrubland Forest or Shrubland Forest or Shrubland
L1: Protected or
IRA Status
Not Protected Not Protected Not Protected
L2: Slope <40% <60% <60%
L3: Distance to
Nearest Road
<1,000ft <2,000ft <2,000ft
L4: Riparian
Buffer
>100ft >100ft >50ft
L5: Administrative
Designation
No Designation No Designation Any Designation

The figure below illustrates the workflow used to quantify the amount of land available for mechanical forest health and risk reduction fuel treatments by considering layered operational constraints. This example illustrates scenario 2 for a 10,252 hectare (25,000 acre) fireshed project area near Rustic, Colorado in the Colorado Front Range priority landscape.

2.2 Data Preparation

2.2.1 Load WCS landscapes

WCS priority landscape spatial data from the USFS Geospatial Data Discovery site (accessed 2023-04-04).

# load data
wf_landscapes <- sf::read_sf("../data/Wildfire_Crisis_Strategy_Landscapes/Wildfire_Crisis_Strategy_Landscapes_(Feature_Layer).shp") %>% 
  dplyr::rename_with(tolower) %>% 
  sf::st_make_valid() %>% 
  dplyr::mutate(
    hectares = (as.numeric(sf::st_area(geometry))/10000) 
    , Mil.Hectares = hectares %>% 
      scales::comma(suffix = " M", scale = 1e-6, accuracy = .01)
    , acres = (as.numeric(sf::st_area(geometry))/4046.85642) 
    , Mil.Acres = acres %>% 
      scales::comma(suffix = " M", scale = 1e-6, accuracy = .01)
  ) %>% 
  dplyr::left_join(
    data.frame(state = datasets::state.name, state_abb = datasets::state.abb)
    , by = join_by("state")
  ) %>% 
  dplyr::mutate(
    area_name = paste0(
      ifelse(is.na(state_abb) | state_abb == "", "UNK", state_abb)
      , ": "
      , name
    )
  ) %>% 
  # manual label placement
  dplyr::left_join(
    readr::read_csv("../data/wildfirepriority/area_label_placement.csv")
    , by = dplyr::join_by("area_name")
  )
#rename sf geom column
  names(wf_landscapes)[names(wf_landscapes)==tolower(attr(wf_landscapes, "sf_column"))] = "geometry"
  sf::st_geometry(wf_landscapes) = "geometry"
# set crs
  transform_crs <- sf::st_crs(wf_landscapes)
# view
  wf_landscapes %>% dplyr::glimpse()

2.2.2 Load landscape-level constraint data

Landscape-level constraint analysis data (i.e., row unique by WCS landscape) was created via this Google Earth Engine script.

# load landscape constraint scenario data
constrained_by_scnro_ls <- list.files("../data/wildfirepriority/all_slp/",pattern = "\\.csv$") %>%
  purrr::keep(stringr::str_starts(.,"wfpriority_all_slp_sc")) %>% 
  purrr::map(function(x){
    readr::read_csv(
      paste0("../data/wildfirepriority/all_slp/",x)
      , name_repair = "universal"
      , col_types = cols(.default = "c")
    ) %>% 
    dplyr::rename_with(tolower) %>% 
    dplyr::rename_with(make.names) %>% 
    dplyr::select(state,name,tidyselect::ends_with("_m2"),tidyselect::starts_with("pct_")) %>% 
    dplyr::mutate(scenario_id = stringr::word(x,4,sep="_") %>% readr::parse_number()) %>% 
    dplyr::relocate(scenario_id)
  }) %>%
  dplyr::bind_rows() %>%
  dplyr::inner_join(
    wf_landscapes %>%
      sf::st_drop_geometry() %>%
      dplyr::select(state,name,area_name)
    , by = dplyr::join_by(state,name)
  ) %>%
  dplyr::relocate(area_name,.after = "scenario_id") %>%
  dplyr::group_by(area_name,scenario_id) %>%
  dplyr::filter(dplyr::row_number()==1) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(
    dplyr::across(
      tidyselect::ends_with("_m2")
      , ~ as.numeric(.x) / 10000
    )
    , dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ as.numeric(.x)
    )
  ) %>%
  dplyr::rename_with(
    ~ gsub("_m2", "_ha", .x)
    , tidyselect::ends_with("_m2")
  ) %>%
  # calculate pct reduction
  dplyr::mutate(
    pct_rdctn1_protected = -1*(1 - pct_rmn1_protected)
    , pct_rdctn2_slope = -1*(pct_rmn1_protected - pct_rmn2_slope)
    , pct_rdctn3_roads = -1*(pct_rmn2_slope - pct_rmn3_roads)
    , pct_rdctn4_riparian = -1*(pct_rmn3_roads - pct_rmn4_riparian)
    , pct_rdctn5_administrative = -1*(pct_rmn4_riparian - pct_rmn5_administrative)
    , pct_rdctn_total = -1*(1 - pct_rmn5_administrative)
    , pct_covertype_area = covertype_area_ha/feature_area_ha
    # scenario description
    , scenario_lab = factor(
      scenario_id
      , levels = 1:3
      , labels = paste0("Scenario ", 1:3)
      , ordered = T
    ) %>% forcats::fct_rev()
    , scenario_desc = factor(
      scenario_id
      , levels = 1:3
      , labels = c("Scenario 1\n(status quo)", "Scenario 2\n(slopes+roads)", "Scenario 3\n(slopes+roads+admin)")
      , ordered = T
    )
  )

2.2.3 Load fireshed project area data

The fireshed registry spatial data was obtained from the USFS Geospatial Data Discovery tool for firesheds and project areas (accessed 2023-05-03).

See this section for exploration of fireshed and project area spatial data

### fireshed spatial data
fireshed <- sf::st_read("../data/firesheds/Fireshed_Registry3A_Fireshed/Fireshed_Registry%3A_Fireshed_(Feature_Layer).shp") %>%
  sf::st_transform(transform_crs) %>% 
  setNames(c(
      "shape_id"
      , "area_ha"
      , "fireshed_id"
      , "fireshed_name"
      , "fireshed_code"
      , "fireshed_state"
      , "nopas"
      , "objectid"
      , "fshed_id"
      , "exp_total"
      , "exp_usfs"
      , "exp_nonfs"
      , "exp_usfs_protected"
      , "exp_nonfs_protected"
      , "exp_usfs_managed"
      , "exp_nonfs_managed"
      , "exp_usfs_forest"
      , "exp_nonfs_forest"
      , "exp_usfs_nonforest"
      , "exp_nonfs_nonforest"
      , "exp_usfs_conifer"
      , "exp_nonfs_conifer"
      , "exp_usfs_managedforest"
      , "exp_nonfs_managedforest"
      , "exp_usfs_managedconifer"
      , "exp_nonfs_managedconifer"
      , "exp_nonfs_nonconifer_hihaz"
      , "dist_vs"
      , "crisis_strategy"
      , "key_preformance_indicator"
      , "national_usfs_rank"
      , "national_all_land_rank"
      , "regional_usfs_rank"
      , "regional_all_land_rank"
      , "start_date"
      , "end_date"
      , "geometry"
  )) %>% 
  dplyr::mutate(
    exposure_pct_rank = dplyr::percent_rank(exp_total)
    , exposure_pct_rank_grp = dplyr::case_when(
      exposure_pct_rank >= 1-0.01 ~ "Top 1%"
      , exposure_pct_rank >= 1-0.05 ~ "Top 5%"
      , exposure_pct_rank >= 1-0.10 ~ "Top 10%"
      , exposure_pct_rank >= 1-0.25 ~ "Top 25%"
      , TRUE ~ "Bottom 75%"
    ) %>% 
    factor(
      levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
      , ordered = T
    )
    # there is also a national_all_land_rank column
    , ntllandrank_pct_rank = dplyr::percent_rank(-national_all_land_rank)
    , ntllandrank_pct_rank_grp = dplyr::case_when(
        ntllandrank_pct_rank >= 1-0.01 ~ "Top 1%"
        , ntllandrank_pct_rank >= 1-0.05 ~ "Top 5%"
        , ntllandrank_pct_rank >= 1-0.10 ~ "Top 10%"
        , ntllandrank_pct_rank >= 1-0.25 ~ "Top 25%"
        , TRUE ~ "Bottom 75%"
      ) %>% 
      factor(
        levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
        , ordered = T
      )
    , crisis_strategy = ifelse(is.na(crisis_strategy),"Not High Risk",crisis_strategy) %>% 
      as.factor() %>% 
      forcats::fct_shift()
  )
  #rename sf geom column
    names(fireshed)[names(fireshed)==tolower(attr(fireshed, "sf_column"))] = "geometry"
    sf::st_geometry(fireshed) = "geometry"
    # calculate area
    fireshed <- fireshed %>% 
      dplyr::mutate(
        fireshed_area_ha = as.numeric(sf::st_area(geometry))/10000
        , fireshed_area_acres = (fireshed_area_ha*10000)/4046.85642
      )
## fireshed_proj_area spatial data
fireshed_proj_area <- sf::st_read("../data/firesheds/Fireshed_Registry3A_Project_Area/Fireshed_Registry%3A_Project_Area_(Feature_Layer).shp") %>%
  sf::st_transform(transform_crs) %>% 
  setNames(c(
      "shape_id"
      , "fireshed_id"
      , "pa_id"
      , "pa_area_ha"
      , "objectid"
      , "pa_id2"
      , "fshed_id"
      , "exp_total"
      , "exp_usfs"
      , "exp_nonfs"
      , "exp_usfs_protected"
      , "exp_nonfs_protected"
      , "exp_usfs_managed"
      , "exp_nonfs_managed"
      , "exp_usfs_forest"
      , "exp_nonfs_forest"
      , "exp_usfs_nonforest"
      , "exp_nonfs_nonforest"
      , "exp_usfs_conifer"
      , "exp_nonfs_conifer"
      , "exp_usfs_managedforest"
      , "exp_nonfs_managedforest"
      , "exp_usfs_managedconifer"
      , "exp_nonfs_managedconifer"
      , "exp_nonfs_nonconifer_hihaz"
      , "dist_vs"
      , "pctrecentlydisturbed"
      , "start_date"
      , "end_date"
      , "geometry"
  )) %>% 
  dplyr::mutate(
    exposure_pct_rank = dplyr::percent_rank(exp_total)
    , exposure_pct_rank_grp = dplyr::case_when(
      exposure_pct_rank >= 1-0.01 ~ "Top 1%"
      , exposure_pct_rank >= 1-0.05 ~ "Top 5%"
      , exposure_pct_rank >= 1-0.10 ~ "Top 10%"
      , exposure_pct_rank >= 1-0.25 ~ "Top 25%"
      , TRUE ~ "Bottom 75%"
    ) %>% 
    factor(
      levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
      , ordered = T
    )
  )
  #rename sf geom column
    names(fireshed_proj_area)[names(fireshed_proj_area)==tolower(attr(fireshed_proj_area, "sf_column"))] = "geometry"
    sf::st_geometry(fireshed_proj_area) = "geometry"
    # calculate area
    fireshed_proj_area <- fireshed_proj_area %>% 
      dplyr::mutate(
        pa_area_ha = as.numeric(sf::st_area(geometry))/10000
        , pa_area_acres = (pa_area_ha*10000)/4046.85642
      ) %>% 
      # JOIN WITH FIRESHED DATA
      dplyr::inner_join(
        fireshed %>%
          sf::st_drop_geometry() %>%
          dplyr::select(fireshed_id, crisis_strategy, exp_total
                        , exposure_pct_rank, exposure_pct_rank_grp
          ) %>% 
          dplyr::rename(exposure_total=exp_total) %>% 
          dplyr::rename_with(
            ~ paste0("fireshed_",.x)
            , -c(fireshed_id)
          )
        , by = dplyr::join_by(fireshed_id)
      ) %>%
      dplyr::select(pa_id,pa_area_ha
                    ,exp_total,exposure_pct_rank,exposure_pct_rank_grp
                    , tidyselect::starts_with("fireshed_")
      ) %>% 
      dplyr::rename(exposure_total=exp_total) %>%
      dplyr::rename_with(
        ~ paste0("pa_", .x, recycle0 = TRUE)
        , tidyselect::starts_with("exp")
      )

Fireshed project area constraint analysis data (i.e., row unique by WCS landscape + f.p.a.) was created via this Google Earth Engine script

# load fireshed constraint scenario data
constrained_by_scnro_ls_pa <-
  list.files("../data/wildfirepriority/fireshed_slp/",pattern = "\\.csv$") %>%
  purrr::keep(stringr::str_starts(.,"wfpriority_fireshed_slp_sc")) %>% 
  purrr::map(function(x){
    readr::read_csv(
      paste0("../data/wildfirepriority/fireshed_slp/",x)
      , name_repair = "universal"
      , col_types = cols(.default = "c")
    ) %>% 
    dplyr::rename_with(tolower) %>% 
    dplyr::rename_with(make.names) %>% 
    dplyr::select(state,name,pa_id,tidyselect::ends_with("_m2"),tidyselect::starts_with("pct_")) %>% 
    dplyr::mutate(scenario_id = stringr::word(x,4,sep="_") %>% readr::parse_number()) %>% 
    dplyr::relocate(scenario_id)
  }) %>%
  dplyr::bind_rows() %>%
  dplyr::inner_join(
    wf_landscapes %>%
      sf::st_drop_geometry() %>%
      dplyr::select(state,name,area_name)
    , by = dplyr::join_by(state,name)
  ) %>%
  dplyr::relocate(area_name,.after = "scenario_id") %>%
  dplyr::mutate(
    dplyr::across(
      tidyselect::ends_with("_m2")
      , ~ as.numeric(.x) / 10000
    )
    , dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ as.numeric(.x)
    )
  ) %>%
  dplyr::rename_with(
    ~ gsub("_m2", "_ha", .x)
    , tidyselect::ends_with("_m2")
  ) %>%
  dplyr::group_by(area_name,pa_id,scenario_id) %>%
  dplyr::arrange(desc(pct_pa_intrsct)) %>%
  dplyr::filter(dplyr::row_number()==1) %>%
  dplyr::ungroup() %>%
  dplyr::filter(pct_pa_intrsct>=0.00001) %>%
  # calculate pct reduction
  dplyr::mutate(
    pct_rdctn1_protected = -1*(1 - pct_rmn1_protected)
    , pct_rdctn2_slope = -1*(pct_rmn1_protected - pct_rmn2_slope)
    , pct_rdctn3_roads = -1*(pct_rmn2_slope - pct_rmn3_roads)
    , pct_rdctn4_riparian = -1*(pct_rmn3_roads - pct_rmn4_riparian)
    , pct_rdctn5_administrative = -1*(pct_rmn4_riparian - pct_rmn5_administrative)
    , pct_rdctn_total = -1*(1 - pct_rmn5_administrative)
    , pct_covertype_area = covertype_area_ha/feature_area_ha
    # calculate level of constraint
      , cnstrnt_lvl = dplyr::case_when(
          -pct_rdctn_total >= 0.85 ~ 1
          , -pct_rdctn_total >= 0.65 ~ 2
          , -pct_rdctn_total >= 0.0 ~ 3
        )
      , cnstrnt_class = factor(
          cnstrnt_lvl 
          , levels = 1:3
          , labels = c("high constraint", "med. constraint", "low constraint")
          , ordered = T
        ) %>% forcats::fct_rev()
      , rmn_cnstrnt_class = factor(
          cnstrnt_lvl 
          , levels = 1:3
          , labels = c("0–15% treatable", "16–35% treatable", ">35% treatable")
          , ordered = T
        ) %>% forcats::fct_rev()
      # scenario description
      , scenario_lab = factor(
          scenario_id
          , levels = 1:3
          , labels = paste0("Scenario ", 1:3)
          , ordered = T
        ) %>% forcats::fct_rev()
      , scenario_desc = factor(
        scenario_id
        , levels = 1:3
        , labels = c("Scenario 1\n(status quo)", "Scenario 2\n(slopes+roads)", "Scenario 3\n(slopes+roads+admin)")
        , ordered = T
      )
  ) %>% 
  # join with fireshed data
  dplyr::inner_join(
    fireshed_proj_area %>%
      sf::st_drop_geometry() %>%
      dplyr::select(pa_id, tidyselect::starts_with("pa_exposure"), tidyselect::starts_with("fireshed_")) %>%
      dplyr::mutate(pa_id=as.character(pa_id))
    , by = dplyr::join_by(pa_id)
  )

2.3 Geographic Delineation Description

Timeline of the research utilized to motivate and support the creation of the WCS:

Nested spatial framework for firesheds. Each scale has specific functionality in terms of the planning processes. Firesheds are the broad scale unit of prioritization, but project areas within them are also prioritized as part of implementation of treatments. Project areas are roughly the size that national forests use for conducting vegetation and fuel management projects. The relative variation among firesheds compared to variation within them controls the relative emphasis on selecting firesheds versus individual planning areas. Ager et al. 2021

#define basemap
states_temp <- USAboundaries::us_states(
  states = c(
        # usfs region 1-6 states
        "MT","WY","CO","NM","AZ","UT","ID","WA","OR","CA","NV"
        # , "KS","NE","SD","ND"
      )
  ) %>% 
  sf::st_transform(transform_crs)
cities_temp <- USAboundaries::us_cities(
  states = c(
        # usfs region 1-6 states
        "MT","WY","CO","NM","AZ","UT","ID","WA","OR","CA","NV"
        # , "KS","NE","SD","ND"
      )
  ) %>% 
  sf::st_transform(transform_crs) %>% 
  dplyr::filter(
    toupper(city) %in% c(
      "DENVER"
      , "TUSCON"
      , "SALT LAKE CITY"
      , "LAS VEGAS"
      , "SAN DIEGO"
      , "LOS ANGELES"
      , "FRESNO"
      , "SAN FRANCISCO"
      , "SACRAMENTO"
      , "PORTLAND"
      , "SEATTLE"
      , "ALBUQUERQUE"
      , "RENO"
      , "BOISE"
      , "HELENA"
      , "BILLINGS"
    )
    | (toupper(city) == "PHOENIX" & state_abbr == "AZ")
  )
# plot basemap
plt_basemap <-
ggplot() + 
  geom_sf(
    data = states_temp
    , fill = NA
    , color = "gray66"
  ) +
  geom_sf_text(
    data = states_temp
    , mapping = aes(label = stusps)
    , color = "gray66"
    , size = 3
  ) +
  geom_sf(
    data = cities_temp
    , shape = 1
    , size = 1
    , color = "gray44"
  ) +
  geom_sf_text(
    data = cities_temp
    , mapping = aes(label = city)
    , color = "gray44"
    , size = 2
    # , hjust = -0.15
    , hjust = 0.5
    , vjust = -0.25
  ) +
  coord_sf(expand = F) +
  theme_void()

# plot basemap simple
plt_basemap_simple <-
ggplot() + 
  geom_sf(
    data = states_temp
    , fill = NA
    , color = "gray66"
  ) +
  coord_sf(expand = F) +
  theme_void()

2.3.1 Fireshed risk to communities (2021)

ggplot() + 
  geom_sf(
    data = fireshed
    , mapping = aes(fill=ntllandrank_pct_rank_grp)
  ) + 
  geom_sf(
    data = USAboundaries::us_states() %>% 
      dplyr::filter(!stusps %in% c("AK","HI","PR")) %>% 
      sf::st_transform(transform_crs)
    , fill = NA
    , color = "gray88"
  ) +
  scale_fill_manual(values=c("firebrick","orange3","khaki","seagreen3","royalblue")) +
  coord_sf(expand = F) +
  labs(
    fill = "Fireshed exposure ranks"
  ) +
  theme_void() + 
  theme(
    legend.position = "top"
    , legend.direction = "horizontal"
    , legend.title = element_text(size = 8)
    , legend.text = element_text(size = 7)
  )

National map of the 7,688 firesheds created from community wildfire transmission data (Evers et al. 2020). The fireshed boundaries were created with a process that delineates hotspots of fire transmission to buildings in adjacent or nearby communities. See the Methods section for details on delineating firesheds. Figure 2 from Ager et al. 2021 (p. 7)

2.3.2 Wilfire Crisis Strategy 21 Priority Landscapes (2022-2023)

plt_fshed <- 
plt_basemap + 
  geom_sf(
    data = fireshed %>%
      sf::st_filter(
        wf_landscapes %>%
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(fill=crisis_strategy)
    , alpha = 0.8
  ) +
  geom_sf(
    data = wf_landscapes
    , mapping = aes(color = paste0(investment, " investment"))
    , fill=NA
    # , color="black" 
    , lwd=0.6
  )+
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = stringr::str_wrap(
          dplyr::case_when(
            stringr::str_count(name, "\\w+")<2
              ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
            , stringr::str_starts(name,"Sierra and Elko") ~ stringr::word(name,1,3)
            , TRUE ~ stringr::word(name,1,2)
          )
          , 7)
        , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , stat = "sf_coordinates"
    , size = 2.5
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
  scale_color_manual(values = c("black", "gray33")) +
  labs(
    fill = "Fireshed\nHigh-Risk Status"
    , color = "Landscapes"
  ) +
  theme(
      legend.position = c(0.87,0.87) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 8)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(lwd = 3, fill = NA)) #, linetype = c(5,1,1)
    , fill = guide_legend(order = 2, override.aes = list(lwd = 0))
  )
# plot
plt_fshed

The boundaries of the 21 priority landscapes roughly follow the boundaries of “firesheds” prioritized to reduce wildfire transmission to developed areas (Figure 1). Firesheds are geographic delineations with an average area of approximately 250,000 acres (~100,000 hectares) that were created to organize the landscape into units for managing wildfire risk to communities.

2.3.3 Wilfire Crisis Strategy Treatment Objective

plt_basemap + 
  geom_sf(
    data = fireshed %>%
      dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
      sf::st_filter(
        wf_landscapes %>%
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(fill=crisis_strategy)
    , alpha = 0.8
  ) +
  geom_sf(
    data = wf_landscapes
    , mapping = aes(color = paste0(investment, " investment"))
    , fill=NA
    # , color="black" 
    , lwd=0.6
  )+
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = stringr::str_wrap(
          dplyr::case_when(
            stringr::str_count(name, "\\w+")<2
              ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
            , stringr::str_starts(name,"Sierra and Elko") ~ stringr::word(name,1,3)
            , TRUE ~ stringr::word(name,1,2)
          )
          , 7)
        , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , stat = "sf_coordinates"
    , size = 2.5
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
  scale_color_manual(values = c("black", "gray33")) +
  labs(
    fill = "Fireshed\nHigh-Risk Status"
    , color = "Landscapes"
  ) +
  theme(
      legend.position = c(0.87,0.87) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 8)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(lwd = 3, fill = NA)) #, linetype = c(5,1,1)
    , fill = guide_legend(order = 2, override.aes = list(lwd = 0))
  )

Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the Strategy (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).

2.3.4 Fireshed Project Areas

name_temp <- "Enchanted Circle"
# pa map
plt_fshed_pa <-
ggplot() + 
  geom_sf(data = 
    fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(
            # fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
            name == name_temp
          ) %>% 
          dplyr::select(pa_id)
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      )
    , aes(fill = fireshed_crisis_strategy, color = "Project Area\nboundary")
    , alpha = 0.8
    , lwd = 1.2
  ) +
  geom_sf(
    data = fireshed %>%
      # dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
      sf::st_filter(
        wf_landscapes %>%
          dplyr::filter(name == name_temp) %>% 
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(color = "Fireshed\nboundary")
    , fill = NA
    , lwd = 1.2
    , linetype = "dashed"
  ) +
  geom_sf(
    data = wf_landscapes %>% dplyr::filter(name == name_temp)
    , mapping = aes(color = "WCS Landscape\nboundary")
    , fill=NA
    , lwd = 0.9
  )+
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
  scale_color_manual(values = c("gray44","skyblue2", "black")) +
  labs(
    fill = "Fireshed\nHigh-Risk Status"
    , color = "Spatial\nFramework"
    , subtitle = name_temp
  ) +
  coord_sf(expand = F) +
  theme_void() +
  theme(
      legend.position = c(1.08,0.79) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 10)
      , plot.subtitle = element_text(face = "bold.italic", size = 10, hjust = 0.5)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(size = 8, linewidth = 5, fill = NA)) #, linetype = c(5,1,1)))
    , fill = guide_legend(order = 2, override.aes = list(size = 8,linewidth = 0, alpha = 1)) # , color = NA
  )
# pa map w/in
plt_fshed_pa_incl <- ggplot() + 
  geom_sf(data = 
    fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(
            # fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
            name == name_temp
          ) %>%
          dplyr::mutate(
            # is_25in = ifelse(pct_pa_intrsct>=0.25,"≥25%","<25%")
            is_25in = ifelse(pct_pa_intrsct>=0.25,"included","excluded")
          ) %>% 
          dplyr::select(
            pa_id
            , is_25in
          )
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      )
    , aes(fill = is_25in, color = "Project Area\nboundary")
    , alpha = 0.9
    , lwd = 1.2
  ) +
  geom_sf(
    data = fireshed %>%
      # dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
      sf::st_filter(
        wf_landscapes %>%
          dplyr::filter(name == name_temp) %>% 
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(color = "Fireshed\nboundary")
    , fill = NA
    , lwd = 1.2
    , linetype = "dashed"
  ) +
  geom_sf(
    data = wf_landscapes %>% dplyr::filter(name == name_temp)
    , mapping = aes(color = "WCS Landscape\nboundary")
    , fill=NA
    , lwd = 0.9
  )+
  scale_fill_manual(values = viridis::viridis(n=4, direction = -1, alpha = 0.9)[c(1,3)]) +
  scale_color_manual(values = c("gray44","skyblue2", "black")) +
  labs(
    fill = "Project Area\nInclusions"
    # fill = "Included in this analysis\nif ≥25% area w/in landscape"
    , color = "Spatial\nFramework"
    , subtitle = name_temp
  ) +
  coord_sf(expand = F) +
  theme_void() +
  theme(
      legend.position = c(1.08,0.83) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 10)
      , plot.subtitle = element_text(face = "bold.italic", size = 10, hjust = 0.5)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(size = 8, linewidth = 5, fill = NA)) #, linetype = c(5,1,1)))
    , fill = guide_legend(order = 2, override.aes = list(size = 8,linewidth = 0, alpha = 1)) # , color = NA
  )

# combine
cowplot::plot_grid(
    plotlist = list(
      NULL
      , plt_fshed_pa + 
        theme(
          plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
          , plot.margin = margin(0.5,3,0.5,0.5, "cm")
        )
      , NULL
      , plt_fshed_pa_incl + 
        theme(
          plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
          , plot.margin = margin(0.5,3,0.5,0.5, "cm")
        )
      , NULL
    )
    , nrow = 1
    , rel_widths = c(0.01,1,0.01,1,0.01)
    , labels = c("","A","","B","")
    , label_size = 17
    , label_fontface = "bold"
    , label_colour = "gray11"
    # , label_y = 0.99
  ) +
  theme(plot.background = element_rect(fill = "white", color = NA))

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_ex_spatial_strctr.png")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 5.8
    , units = "in"
    , dpi = "print"
  )
  1. Nested within firesheds are project areas of approximately 10,000 ha in size which represent the geographic unit at which vegetation and fuel management projects are planned (Ager et al. 2021).

  2. Fireshed project areas that did not have at least 25% of area within the boundary of the Strategy priority landscape were excluded from this analysis. This cutoff was implemented under the assumption that with less than 25% of area available, treatment alone would not substantially affect wildfire behavior across the fireshed (North et al. 2015). For each fireshed project area included in this analysis, the entire area within the project area boundary was used to calculate the percent of mechanically available acreage including area outside of the priority landscape boundary.

2.4 Landscape-level analysis

# filter fireshed project area
constrained_by_scnro_ls_pa <- constrained_by_scnro_ls_pa %>% 
  dplyr::filter(pct_pa_intrsct>=0.25)

2.4.1 Reduction Treatable Area Table

# scnum_temp<-2
  tbl_temp <- constrained_by_scnro_ls %>% 
    # dplyr::filter(scenario_id==scnum_temp) %>% 
    dplyr::mutate(
      dplyr::across(
        tidyselect::ends_with("_ha")
        , ~ scales::comma(.x, accuracy = 1)
      )
      , dplyr::across(
        tidyselect::starts_with("pct_")
        , ~ scales::percent(.x, accuracy = 0.1)
      )
    ) %>% 
    dplyr::select(
      area_name
      , scenario_lab
      , covertype_area_ha
      , pct_covertype_area
      , pct_rdctn1_protected
      , pct_rdctn2_slope
      , pct_rdctn3_roads
      , pct_rdctn4_riparian
      , pct_rdctn5_administrative
      , rmn5_administrative_area_ha
      , pct_rmn5_administrative
    ) %>% 
    dplyr::arrange(area_name,desc(scenario_lab))
  # make table 
  kableExtra::kable(
      tbl_temp %>% dplyr::select(-c(area_name))
      # , caption = paste0("<b><font color=navy>Scenario "
      #                    ,scnum_temp
      #                    ,"</font></b><br>percent reduction of different types of constraints on mechanical treatment"
      #         )
      , caption = "<b>Wildfire Crisis Strategy 21 Priority Landscapes</b><br>percent reduction of different types of constraints on mechanical treatment"
      , col.names = c(
        ""
        , "Forest+Shrub\n(ha)", "Forest+Shrub\n%"
        , "Protected"
        , "Slope\nSteepness"
        , "Road\nDistance"
        , "Riparian\nBuffer"
        , "Administrative"
        , "Remaining (ha)"
        , "Remaining (%)"
      )
      , escape = F
    ) %>%  
    add_header_above(c(" " = 1, "Area of\nPriority Landscape"=2, "Constraint\nLeast Flexible to Most Flexible" = 5, " " = 2)) %>%
    kable_classic(full_width=T) %>% 
    pack_rows(index = table(forcats::fct_inorder(tbl_temp$area_name))) %>% 
    kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%  
    kableExtra::scroll_box(width = "740px")
Table 2.2: Wildfire Crisis Strategy 21 Priority Landscapes
percent reduction of different types of constraints on mechanical treatment
Area of
Priority Landscape
Constraint
Least Flexible to Most Flexible
Forest+Shrub (ha) Forest+Shrub % Protected Slope Steepness Road Distance Riparian Buffer Administrative Remaining (ha) Remaining (%)
AZ: 4FRI
Scenario 1 2,190,804 90.6% -13.8% -7.5% -26.4% -10.0% -15.2% 592,541 27.0%
Scenario 2 2,190,804 90.6% -13.8% -2.0% -13.2% -13.9% -19.7% 819,126 37.4%
Scenario 3 2,190,804 90.6% -13.8% -2.0% -13.2% -9.8% 0.0% 1,339,032 61.1%
AZ: Prescott
Scenario 1 219,110 86.4% -12.1% -16.6% -24.0% -12.1% -5.3% 65,793 30.0%
Scenario 2 219,110 86.4% -12.1% -3.2% -12.2% -18.9% -8.3% 99,404 45.4%
Scenario 3 219,110 86.4% -12.1% -3.2% -12.2% -13.4% 0.0% 129,550 59.1%
AZ: San Carlos Apache Tribal Forest Protection
Scenario 1 1,056,065 91.4% -6.1% -20.4% -51.9% -5.8% -0.8% 157,750 14.9%
Scenario 2 1,056,065 91.4% -6.1% -5.3% -48.8% -10.2% -1.7% 294,497 27.9%
Scenario 3 1,056,065 91.4% -6.1% -5.3% -48.8% -7.3% 0.0% 342,834 32.5%
CA: North Yuba
Scenario 1 139,202 96.5% -11.6% -33.4% -7.8% -14.7% -2.3% 42,013 30.2%
Scenario 2 139,202 96.5% -11.6% -13.9% -3.2% -23.0% -3.3% 62,696 45.0%
Scenario 3 139,202 96.5% -11.6% -13.9% -3.2% -16.7% 0.0% 76,029 54.6%
CA: Plumas Community Protection
Scenario 1 100,872 93.3% 0.0% -23.2% -19.8% -21.5% -1.9% 33,888 33.6%
Scenario 2 100,872 93.3% 0.0% -5.2% -8.1% -33.9% -2.9% 50,328 49.9%
Scenario 3 100,872 93.3% 0.0% -5.2% -8.1% -24.4% 0.0% 62,810 62.3%
CA: Southern California Fireshed Risk Reduction Strategy
Scenario 1 1,408,168 86.0% -57.0% -21.0% -8.3% -2.3% -1.4% 139,165 9.9%
Scenario 2 1,408,168 86.0% -57.0% -9.6% -6.5% -4.2% -3.0% 277,540 19.7%
Scenario 3 1,408,168 86.0% -57.0% -9.6% -6.5% -3.0% 0.0% 336,002 23.9%
CA: Stanislaus
Scenario 1 116,139 94.3% -2.9% -26.7% -10.6% -24.7% -0.6% 40,057 34.5%
Scenario 2 116,139 94.3% -2.9% -10.7% -3.9% -35.8% -0.8% 53,380 46.0%
Scenario 3 116,139 94.3% -2.9% -10.7% -3.9% -26.1% 0.0% 65,599 56.5%
CA: Trinity Forest Health and Fire-Resilient Rural Communities
Scenario 1 295,916 80.4% -12.9% -39.2% -10.8% -10.7% -12.0% 42,598 14.4%
Scenario 2 295,916 80.4% -12.9% -11.9% -6.5% -20.7% -24.1% 70,790 23.9%
Scenario 3 295,916 80.4% -12.9% -11.9% -6.5% -15.1% 0.0% 158,737 53.6%
CO: Colorado Front Range
Scenario 1 1,051,540 72.6% -20.4% -19.9% -20.0% -11.3% -1.6% 281,048 26.7%
Scenario 2 1,051,540 72.6% -20.4% -5.6% -12.0% -17.3% -3.2% 435,578 41.4%
Scenario 3 1,051,540 72.6% -20.4% -5.6% -12.0% -12.3% 0.0% 521,905 49.6%
ID: Nez Perce-Clearwater-Lower Salmon
Scenario 1 699,884 89.0% -43.0% -19.7% -9.1% -1.4% -0.3% 185,161 26.5%
Scenario 2 699,884 89.0% -43.0% -6.7% -5.5% -2.3% -0.7% 292,841 41.8%
Scenario 3 699,884 89.0% -43.0% -6.7% -5.5% -1.6% 0.0% 302,267 43.2%
ID: Southwest Idaho
Scenario 1 608,338 87.6% -13.3% -28.3% -10.4% -2.0% -0.3% 278,231 45.7%
Scenario 2 608,338 87.6% -13.3% -7.0% -6.6% -2.8% -0.8% 422,974 69.5%
Scenario 3 608,338 87.6% -13.3% -7.0% -6.6% -2.0% 0.0% 432,970 71.2%
MT: Kootenai Complex
Scenario 1 372,345 91.4% -21.4% -15.0% -7.8% -6.9% -10.6% 142,439 38.3%
Scenario 2 372,345 91.4% -21.4% -3.6% -2.1% -8.6% -13.8% 188,288 50.6%
Scenario 3 372,345 91.4% -21.4% -3.6% -2.1% -6.1% 0.0% 249,151 66.9%
NM: Enchanted Circle
Scenario 1 506,460 85.6% -8.5% -21.5% -33.4% -7.5% -2.2% 135,772 26.8%
Scenario 2 506,460 85.6% -8.5% -5.2% -25.7% -11.7% -3.7% 228,840 45.2%
Scenario 3 506,460 85.6% -8.5% -5.2% -25.7% -8.2% 0.0% 265,124 52.3%
NV: Sierra and Elko Fronts
Scenario 1 1,003,525 73.2% -26.6% -9.8% -26.9% -6.2% -0.9% 297,023 29.6%
Scenario 2 1,003,525 73.2% -26.6% -2.0% -15.8% -8.9% -1.5% 454,449 45.3%
Scenario 3 1,003,525 73.2% -26.6% -2.0% -15.8% -6.3% 0.0% 495,822 49.4%
OR: Central Oregon
Scenario 1 937,977 88.8% -7.2% -2.7% -14.9% -5.1% -10.2% 562,234 59.9%
Scenario 2 937,977 88.8% -7.2% -0.5% -4.7% -6.0% -11.6% 656,297 70.0%
Scenario 3 937,977 88.8% -7.2% -0.5% -4.7% -4.2% 0.0% 781,760 83.3%
OR: Klamath River Basin
Scenario 1 2,900,530 77.7% -20.5% -16.9% -17.0% -5.5% -4.6% 1,032,327 35.6%
Scenario 2 2,900,530 77.7% -20.5% -5.7% -8.1% -8.6% -8.2% 1,420,954 49.0%
Scenario 3 2,900,530 77.7% -20.5% -5.7% -8.1% -6.2% 0.0% 1,729,729 59.6%
OR: Mount Hood Forest Health and Fire-Resilient Communities
Scenario 1 334,397 76.7% -19.6% -14.2% -16.9% -9.1% -16.4% 79,550 23.8%
Scenario 2 334,397 76.7% -19.6% -4.1% -7.3% -13.2% -20.9% 116,660 34.9%
Scenario 3 334,397 76.7% -19.6% -4.1% -7.3% -9.5% 0.0% 198,932 59.5%
UT: Pine Valley
Scenario 1 153,246 94.2% -37.2% -6.9% -21.5% -6.7% 0.0% 42,532 27.8%
Scenario 2 153,246 94.2% -37.2% -1.4% -10.5% -9.2% 0.0% 63,919 41.7%
Scenario 3 153,246 94.2% -37.2% -1.4% -10.5% -6.5% 0.0% 68,040 44.4%
UT: Wasatch
Scenario 1 381,429 89.5% -43.0% -12.1% -14.8% -4.9% -0.3% 95,232 25.0%
Scenario 2 381,429 89.5% -43.0% -3.7% -7.9% -6.7% -0.5% 145,982 38.3%
Scenario 3 381,429 89.5% -43.0% -3.7% -7.9% -4.8% 0.0% 155,181 40.7%
WA: Central Washington Initiative
Scenario 1 670,440 68.1% -20.6% -32.0% -12.3% -8.5% -12.2% 96,766 14.4%
Scenario 2 670,440 68.1% -20.6% -10.4% -7.8% -14.1% -21.8% 168,986 25.2%
Scenario 3 670,440 68.1% -20.6% -10.4% -7.8% -10.2% 0.0% 341,540 50.9%
WA: Colville Northeast Washington Vision
Scenario 1 655,989 92.8% -20.1% -18.5% -16.0% -10.1% -1.2% 223,368 34.1%
Scenario 2 655,989 92.8% -20.1% -3.8% -7.5% -14.4% -2.0% 342,407 52.2%
Scenario 3 655,989 92.8% -20.1% -3.8% -7.5% -10.4% 0.0% 382,335 58.3%

2.4.2 Change in Available by Scenario

constrained_by_scnro_ls %>% 
  dplyr::mutate(scenario_lab = scenario_lab %>% forcats::fct_rev()) %>% 
ggplot(
    mapping = aes(x=scenario_lab,y=pct_rmn5_administrative,group=area_name)
    ) +
    geom_line(mapping=aes(color = area_name), linewidth = 1.5) +
    geom_label(
      mapping=aes(label = scales::percent(pct_rmn5_administrative, scale = 100, accuracy = 1))
      , color = "black"
      , size = 3
      , label.padding = unit(0.15, "lines")
    ) +
    facet_wrap(facets = vars(area_name)
      , ncol = 3
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
      , scales = "free_x"
    ) +
    scale_color_viridis_d(option = "turbo", alpha = 0.8) +
    scale_y_continuous(expand = expansion(mult = c(0.1, 0.1)), labels = scales::percent_format(), breaks = scales::extended_breaks(6)) +
    # scale_x_discrete(labels = scales::label_wrap(20)) +
    labs(
      x = "Most Restrictive \U2192 Least Restrictive"
      , y = "Percent Treatable Area Remaining"
    ) +
    theme_light() +
    theme(
      legend.position = "none"
      , axis.text.y = element_text(size=7)
      , axis.text.x = element_text(size=7.5)
      , strip.text = element_text(color = "black", face = "bold", size = 10)
    )

2.4.3 Reduction by Constraint

# reshape
constrained_by_scnro_ls_long <- constrained_by_scnro_ls %>%
  dplyr::select(state, name, area_name, tidyselect::starts_with("scenario_"), tidyselect::starts_with("pct_rdctn")) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("pct_rdctn")
    , names_to = "constraint"
    , values_to = "pct_rdctn"
    , names_prefix = "pct_rdctn"
    , values_drop_na = F
  ) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , constraint = factor(
        constraint
        , ordered = TRUE
        , levels = c(
          "protected"
          , "slope"
          , "roads"
          , "riparian"
          , "administrative"
          , "total"   
        )
        , labels = c(
          "Protected"
          , "Slope\nSteepness"
          , "Road\nDistance"
          , "Riparian\nBuffer"
          , "Administrative"
          , "Total"
        )
      ) %>% forcats::fct_rev()
  ) %>% 
  dplyr::left_join(
    constrained_by_scnro_ls %>% dplyr::select(state,name,scenario_id,pct_rdctn_total)
    , by = join_by(state,name,scenario_id)
  )
# color pallete removing dark blue
n_temp <- ((constrained_by_scnro_ls_long$constraint %>% unique() %>% length())-1)*2
plasma_custom <- viridis::plasma(n = n_temp)[seq(2,n_temp,2)]
# scales::show_col(plasma_custom)
# plasma_custom
# plot
ggplot(
  data = constrained_by_scnro_ls_long %>% 
      dplyr::filter(constraint!="Total") %>% 
      dplyr::mutate(
        dplyr::across(
          c(scenario_lab, area_name)
          , ~ forcats::fct_rev(.x)
        )
      )
  , mapping = aes(y = area_name, x = pct_rdctn)
) +
  geom_col(
    mapping = aes(fill = constraint)
    , color = NA, width = 0.7
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      label = scales::percent(ifelse(pct_rdctn<0.11*-1,pct_rdctn,NA), accuracy = 1)
      , group = constraint
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_long %>% 
      dplyr::filter(constraint=="Total") %>% 
      dplyr::mutate(
        dplyr::across(
          c(scenario_lab, area_name)
          , ~ forcats::fct_rev(.x)
        )
      )
    , mapping = aes(
      y = area_name, x = pct_rdctn_total
      , label = scales::percent(pct_rdctn_total, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_grid(cols = vars(scenario_lab)) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0.04, .3)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x =  element_blank() #
    , strip.text = element_text(color = "black", face = "bold")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

ggplot2::ggsave(
    filename = paste0("../data/plt_cnstrnt_rdctn_ls.png")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

2.5 Fireshed-level analysis

Based on model simulations of how much area generally needs to be treated to influence wildfire behavior, we binned the firesheds into three classes of mechanical constraint:

  • high (85–100% [i.e., only 0–15% is available for mechanical treatment]): fuels treatment would principally need to rely on fire
  • medium (65– 84%): could use a combination of fire and mechanical thinning
  • low (65%): could effectively influence wildfire behavior with mechanical treatment alone

2.5.1 Distribution of Fireshed PA Constraint

constrained_by_scnro_ls_pa %>% 
  dplyr::count(area_name,scenario_lab,cnstrnt_class) %>% 
  dplyr::group_by(area_name,scenario_lab) %>% 
  dplyr::mutate(
    pct = n/sum(n)
    , high_pct=max(ifelse(cnstrnt_class=="high constraint",pct,0)) 
    , tot = sum(n)
    , scenario_lab = scenario_lab %>% forcats::fct_rev()
  ) %>% 
ggplot() +
  geom_col(
    mapping = aes(x = pct, y = reorder(area_name, desc(area_name)), fill=cnstrnt_class)
    ,width = 0.7, alpha=0.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = reorder(area_name, desc(area_name)), group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.055,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  # geom_text(
  #     data = constrained_byftr_huc12_wide %>%
  #       dplyr::group_by(area_name) %>%
  #       dplyr::summarise(n=n(), high_n=sum(ifelse(cnstrnt_class=="high constraint",1,0)))
  #     , mapping = aes(x = -0.05,y=reorder(area_name, -(high_n/n)),label=paste0("n=",scales::comma(n)))
  #     , size = 3
  #     , color = "black"
  #     , hjust = 0.7
  #   ) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  facet_grid(cols = vars(scenario_lab)) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "% of Fireshed\nProject Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

2.5.2 Distribution of Fireshed Area

# plot
ggplot(data = constrained_by_scnro_ls_pa %>% dplyr::filter(scenario_id==1)
         , mapping = aes(x = feature_area_ha, y = reorder(area_name,desc(area_name) ))
    ) +
    geom_vline(
      xintercept = 
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(scenario_id==1) %>% 
          dplyr::pull(feature_area_ha) %>% 
          median()
      , linetype="dashed", color="gray66"
    ) +
    geom_violin(aes(fill = state)) + 
    geom_boxplot(width = 0.15) +
    # geom_point(size = 0.7, color = "gray33", alpha = 0.2) +
    geom_text(
      data = constrained_by_scnro_ls_pa %>% 
        dplyr::filter(scenario_id==1) %>% 
        dplyr::group_by(area_name) %>% 
        dplyr::summarise(n=n(), max_area=max(feature_area_ha))
      , mapping = aes(
          x = (
            constrained_by_scnro_ls_pa %>% 
            dplyr::filter(scenario_id==1) %>% 
            dplyr::pull(feature_area_ha) %>% 
            min()
          )*.9
          ,y=reorder(area_name,desc(area_name) )
          ,label=paste0("n=",scales::comma(n))
        )
      , size = 3
      , color = "black"
      , hjust = 0.7
    ) +
    scale_fill_viridis_d(option = "viridis", alpha = 0.8) +
    scale_x_continuous(labels = scales::comma, breaks = scales::extended_breaks(n=7)) +
    labs(
      fill = ""
      , y = ""
      , x = "Fireshed Area (ha)"
    ) +
    theme_light() +
    theme(
      legend.position = "none" # c(0.9, 0.9)
      , axis.title = element_text(size=9)
      , axis.text.x = element_text(size=7)
    )

2.5.3 Map of Firesheds by Level of Constraint

# polish for mapping
map_fs_temp <- fireshed_proj_area %>% 
  dplyr::select(pa_id, geometry) %>% 
  dplyr::mutate(pa_id=as.character(pa_id)) %>% 
  dplyr::inner_join(
    constrained_by_scnro_ls_pa
    , by = dplyr::join_by("pa_id")
    , multiple = "all"
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ scales::percent(as.numeric(.x), accuracy = 0.1)
    )
    , dplyr::across(
      tidyselect::ends_with("_ha")
      , ~ scales::comma(as.numeric(.x), accuracy = 1)
    )
    , scenario_lab = scenario_lab %>% forcats::fct_rev()
  ) %>% 
  dplyr::mutate(
    Priority.Area.Name = area_name
    , Level.of.Constraint = cnstrnt_class
    , Pct.Treatable.Class = rmn_cnstrnt_class
    , Fireshed.ProjArea.ID = pa_id
    , Fireshed.ProjArea.Area.Hectares = feature_area_ha
    , Fireshed.WCS.Status = ifelse(is.na(fireshed_crisis_strategy),"Not High Risk",fireshed_crisis_strategy)
    , ForestShrub.Area.Hectares = covertype_area_ha
    , Pct.Rdctn.Protected = pct_rdctn1_protected
    , Pct.Rdctn.Slope = pct_rdctn2_slope
    , Pct.Rdctn.Roads = pct_rdctn3_roads
    , Pct.Rdctn.Riparian = pct_rdctn4_riparian
    , Pct.Rdctn.Administrative = pct_rdctn5_administrative
    , Treatable.Area.Hectares = rmn5_administrative_area_ha
    , Pct.Treatable.Area = pct_rmn5_administrative
  )
# filter for mapview
mapview_fs_temp <- map_fs_temp %>% 
  dplyr::filter(
    scenario_id == 2
  )
# ## filter for high priority only
# hi_map_fs_temp <- map_fs_temp %>% 
#   dplyr::filter(!is.na(fireshed_crisis_strategy))
# filter for mapview
mapview_fs_sc1_temp <- map_fs_temp %>% 
  dplyr::filter(
    scenario_id == 1
  )
mapview_fs_sc3_temp <- map_fs_temp %>% 
  dplyr::filter(
    scenario_id == 3
  )
mapview_fs_sc3_hi_temp <- map_fs_temp %>% 
  dplyr::filter(
    scenario_id == 3
    & fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
  )
# mapview
mapview::mapviewOptions(homebutton = FALSE, basemaps = c("OpenStreetMap","Esri.WorldImagery"))
mapview::mapview(wf_landscapes
        , color = "black"
        , lwd = 1
        , alpha.regions = 0
        , label = FALSE
        , legend = FALSE
        , popup = FALSE
        , layer.name = "WCS Landscapes"
        , hide = TRUE
  ) +
mapview::mapview(
  mapview_fs_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 2 Constraints"
  , hide = FALSE
  , popup = leafpop::popupTable(
      mapview_fs_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) + 
mapview::mapview(
  mapview_fs_sc1_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 1 Constraints"
  , hide = TRUE
  , popup = leafpop::popupTable(
      mapview_fs_sc1_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) + 
mapview::mapview(
  mapview_fs_sc3_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 3 Constraints"
  , hide = TRUE
  , popup = leafpop::popupTable(
      mapview_fs_sc3_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) +
mapview::mapview(
  mapview_fs_sc3_hi_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 3 Constraints (High Only)"
  , hide = TRUE
  , popup = leafpop::popupTable(
      mapview_fs_sc3_hi_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) 
# **Scenario 2 used for map above**

2.5.4 Map of Fireshed Project Areas by Scenario

plt_fn_temp <- function(scnum=1) {
  plt_basemap_simple +
    geom_sf(data = map_fs_temp %>% dplyr::filter(scenario_id==scnum)
            , aes(fill = cnstrnt_class), lwd = 0.05) +
    
  geom_sf(
    data = wf_landscapes
    # , mapping = aes(color = paste0(investment, " investment"))
    , fill=NA
    , color="black"
    , lwd=0.2
  )+
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = stringr::str_wrap(
          dplyr::case_when(
            stringr::str_count(name, "\\w+")<2
              ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
            , stringr::str_starts(name,"Sierra and Elko") ~ stringr::word(name,1,3)
            , TRUE ~ stringr::word(name,1,2)
          )
          , 7)
        # , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , color = "black"
    , stat = "sf_coordinates"
    , size = 2
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  # scale_color_manual(values = c("black", "gray33")) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  labs(
    fill = ""
    , tag = paste0(
      map_fs_temp %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::pull(scenario_lab) %>% unique()
      , " Map"
    )
  ) +
  theme(
    legend.position = c(0.8,0.9) # "top"
    , legend.text = element_text(size = 7)
    , plot.title = element_text(face = "bold", size = 9, margin = margin(0,0,30,0))
    , plot.tag = element_text(face = "bold", size = 10)
    , plot.tag.position = c(0.1,0.98) # "top"
    # , plot.tag.position = c(0.8,0.94) # "top"
  ) + 
  guides(color = "none")
}
# plot
# plt_fn_temp()
plt_fpas_sc_map <- 
  map_fs_temp %>% dplyr::pull(scenario_id) %>% unique() %>% 
    purrr::map(plt_fn_temp)
plt_fpas_sc_map

2.5.5 Map Fireshed Project Areas by Scenario and Landscape

How do small multiple maps of scenarios look?

the answer is: bad see this post about small multiple maps

area_nm_list <- sort(unique(constrained_by_scnro_ls$area_name))
facet_map_fn <- function(row_n=1) {
  plt <- ggplot() +
    geom_sf(
      data = map_fs_temp %>% 
        dplyr::filter(
          area_name == area_nm_list[row_n]
        )
      , aes(fill = cnstrnt_class), lwd = 0.05
    ) +
    geom_sf(
      data = wf_landscapes %>% 
        dplyr::filter(
          area_name == area_nm_list[row_n]
        )
      , color = "black"
      , fill=NA
      , lwd=0.2
    ) +
    facet_grid(
      cols = vars(scenario_lab)
      , rows = vars(area_name)
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
      , switch = "both"
      # , switch = "x"
    ) +
    scale_color_manual(values = c("black", "gray33")) +
    scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
    labs(
      fill = ""
      , tag = stringr::str_wrap(area_nm_list[row_n], 35)
        # scales::label_wrap(area_nm_list[row_n], 10)
    ) +
    theme_void() +
    # theme_light() +
    theme(
      legend.position = "none"
      # , legend.text = element_text(size = 7)
      , plot.subtitle = element_text(face = "bold", size = 9, hjust = 0.5, vjust = -5)
      # attempting to fix facet labels not aligning
      , strip.text.x = element_text(size = 8, vjust = 1)
      , strip.text.y = element_blank()
      , plot.tag.position = c(.5, 0.98)
      , plot.tag = element_text(size = 8, face = "bold")
      # , plot.background = element_rect(color = "gray88", fill=NA, size=0.5)
      , plot.margin = margin(0, 0, 0, 0, "cm")
      # facet labels not alinging in plot_grid
      # , strip.text.x = element_text(color = "black", face = "bold")
      # , strip.text.y.left = element_text(angle = 0)
      # , strip.placement = "outside"
      
    ) + 
    guides(color = "none")
  return(plt)
}
# print(facet_map_fn())
plt_list <- 1:length(area_nm_list) %>% 
  purrr::map(facet_map_fn)
# combine
scenario_maps <- cowplot::plot_grid(
  plotlist = plt_list #[1:4]
  , ncol = 3
  , align = "hv"
)
# scenario_maps
# plot_grid didn't work ;(
plt_list[16]
# plt_list

Make each plot individually by landscape

#########################################################################
#########################################################################
#########################################################################
##################hack to align plots for ggmap
ggmap_bbox_fn <- function(map, my_crs=3857) {
    if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
    # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
    # and set the names to what sf::st_bbox expects:
    map_bbox <- setNames(unlist(attr(map, "bb")), c("ymin", "xmin", "ymax", "xmax"))
    # Convert the bbox to an sf polygon, transform it to 3857, 
    # and convert back to a bbox (convoluted, but it works)
    bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), my_crs))
    # Overwrite the bbox of the ggmap object with the transformed coordinates 
    attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
    attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
    attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
    attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
    map
}
plt_crs <- 3857
#########################################################################
#########################################################################
#########################################################################
area_nm_list <- sort(unique(constrained_by_scnro_ls$area_name))
# buffer for smaller areas
if_smaller_ha_temp <- wf_landscapes %>% 
  dplyr::pull(hectares) %>% 
  quantile(0.45) %>% 
  purrr::pluck(1)
# create color palette
cols_rdylbu_lt <- c(
  "low constraint" = rev(RColorBrewer::brewer.pal(n=3,name="RdYlBu"))[1]
  , "med. constraint" = rev(RColorBrewer::brewer.pal(n=3,name="RdYlBu"))[2]
  , "high constraint" = rev(RColorBrewer::brewer.pal(n=3,name="RdYlBu"))[3]
)
# function to plot basemap
landscape_basemap_fn <- function(row_n=1) {
  # should zoom in?
  zoom_level <- dplyr::case_when(
    area_nm_list[row_n] %in% c(
      "WA: Colville Northeast Washington Vision"
      , "CA: Southern California Fireshed Risk Reduction Strategy"
      ) ~ 7
    , area_nm_list[row_n] %in% c(
      "UT: Wasatch"
      , "CA: Trinity Forest Health and Fire-Resilient Rural Communities"
      , "ID: Nez Perce-Clearwater-Lower Salmon"
      , "NV: Sierra and Elko Fronts"
      , "OR: Mount Hood Forest Health and Fire-Resilient Communities"
    ) ~ 8
    , area_nm_list[row_n] %in% c(
      "CA: Plumas Community Protection"
      , "CA: North Yuba"
      , "MT: Kootenai Complex"
      , "UT: Pine Valley"
    ) ~ 9
    , 
      (wf_landscapes %>% 
        dplyr::filter(
            area_name == area_nm_list[row_n]
        ) %>% 
        dplyr::pull(hectares) %>% 
        purrr::pluck(1)
        <= if_smaller_ha_temp
      ) ~ 10
    , TRUE ~ 8
  )
  # should buffer extend?
  buffer_box <- dplyr::case_when(
    area_nm_list[row_n] %in% c("WA: Colville Northeast Washington Vision") ~ 50000
    , area_nm_list[row_n] %in% c(
      "CA: Southern California Fireshed Risk Reduction Strategy"
      , "ID: Nez Perce-Clearwater-Lower Salmon"
      , "MT: Kootenai Complex"
    ) ~ 40000
    , area_nm_list[row_n] %in% c(
      "AZ: San Carlos Apache Tribal Forest Protection"
      , "NV: Sierra and Elko Fronts"
      , "UT: Pine Valley"
    ) ~ 21000
    , area_nm_list[row_n] %in% c(
      "CA: Plumas Community Protection"
      , "CA: North Yuba"
      , "CA: Trinity Forest Health and Fire-Resilient Rural Communities"
      , "OR: Mount Hood Forest Health and Fire-Resilient Communities"
    ) ~ 35000
    , area_nm_list[row_n] %in% c("UT: Wasatch") ~ 5000
    , 
      (wf_landscapes %>% 
        dplyr::filter(
            area_name == area_nm_list[row_n]
        ) %>% 
        dplyr::pull(hectares) %>% 
        purrr::pluck(1)
        <= if_smaller_ha_temp
      ) ~ 15000
    , TRUE ~ 5000
  )
  
  # bounding box
  bb_temp <- 
    map_fs_temp %>% 
    dplyr::filter(
      area_name == area_nm_list[row_n]
    ) %>% 
    sf::st_union() %>% 
    sf::st_transform(crs=5070) %>% 
    sf::st_buffer(as.numeric(buffer_box)) %>% 
    sf::st_transform(crs=4326) %>% # same as get_map return
    sf::st_bbox()
  # set bbox for get call
  bbox_temp <- c(
    bottom = bb_temp[[2]]
    , top = bb_temp[[4]]
    , right = bb_temp[[3]]
    , left = bb_temp[[1]]
  )
  # get map
  hey_ggmap <- ggmap::get_stamenmap(
    bbox = bbox_temp
    , zoom = zoom_level
    , maptype = "terrain" #"toner-hybrid" # "terrain"
    , crop = T
  )
  # ggmap(hey_ggmap)
  # apply align function
    hey_ggmap_aligned <- ggmap_bbox_fn(hey_ggmap, plt_crs) # Use the function
  # plot
  plt <- ggmap(hey_ggmap_aligned) + 
    coord_sf(
      expand = FALSE
    ) +
    labs(
      title = area_nm_list[row_n]
    ) +
    theme_light() +
    theme(
      legend.position = "bottom"
      , legend.direction  = "horizontal"
      , legend.title = element_text(size=7)
      , legend.margin = margin(c(-10,0,0,0))
      , plot.title = element_text(size = 11)
      , strip.text = element_text(color = "black", face = "bold")
      , axis.title = element_blank()
      , axis.text = element_blank()
      , axis.ticks = element_blank()
      , panel.grid = element_blank()
      , plot.margin = margin(0, 0, 0, 0, "cm")
    )
  # return
  return(plt)
}
# landscape_basemap_fn(row_n = 9)
# function to plot scenario constraints + basemap
scenario_cnstrnt_map_fn <- function(row_n = 1) {
  plt <- landscape_basemap_fn(row_n) + 
      geom_sf(
        data = map_fs_temp %>%
          dplyr::filter(
            area_name == area_nm_list[row_n]
          ) %>%
          sf::st_transform(crs=plt_crs)
        , aes(fill = cnstrnt_class), lwd = 0.05
        , inherit.aes = F
        , alpha = 0.6
      ) +
      geom_sf(
        data = wf_landscapes %>% 
          dplyr::filter(
              area_name == area_nm_list[row_n]
            ) %>% 
          sf::st_transform(crs=plt_crs)
        , fill = NA, color = "black", lwd = 0.3
        , inherit.aes = F
      ) +
      facet_grid(
        cols = vars(scenario_lab)
        , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
        # , switch = "both"
      ) +
      scale_fill_manual(values = cols_rdylbu_lt) +
      labs(
        fill = "Level of\nConstraint"
      ) +
        guides(
        fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
      )
  # return
  return(plt)
}
# scenario_cnstrnt_map_fn(row_n = 9)
# go through list of areas to make plot
plt_list <- 1:length(area_nm_list) %>%
  purrr::map(scenario_cnstrnt_map_fn)
# scenario_maps
plt_list[]

####
# area_nm_list
# scenario_cnstrnt_map_fn(18)

2.5.6 Patch Analysis

Notes from Collins et al. 2010:

The larger the individual fuel treatment the greater the potential for tree survival, because of reduced edge effect mortality (Ritchie et al. 2007). Larger individual treatments have a greater potential to reduce fire behavior and slow fire spread, which ultimately impacts adjacent untreated stands and should enhance suppression opportunities and increase firefighter safety (Collins et al. 2010). Planning fewer, larger individual treatments across the landscape appears to be a better strategy when human community protection is a primary concern (Schmidt et al. 2008). Finney et al. (2007) note that the relative improvement of optimized treatment placement breaks down when larger proportions of the landscape (40–50%) are excluded from treatment because of land-management constraints.

# create spatial data of interconnected patches by constraint class
# this is to represent the spatial aggregation of fireshed project areas of similar levels of  mechanical constraint
  # aggregation = Tendency of patch or land-cover types to be spatially adjacent or in close proximity
cnstrnt_class_patches <-
  fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::select(pa_id, area_name, scenario_id, scenario_lab, cnstrnt_class)
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      ) %>% 
  # union vectors/polygons by constraint class into one multipolygon
    # this dissolves inner boundaries by constraint class
  dplyr::group_by(area_name, scenario_id, scenario_lab, cnstrnt_class) %>% 
  dplyr::summarize(
    geometry = sf::st_union(geometry, by_feature = F)
  ) %>% 
  # separate a multipolygon geometry into several polygons objects after performing a st_union()
  sf::st_cast(to = "MULTIPOLYGON") %>%
  st_cast(to = "POLYGON") %>% 
  # calculate patch area 
  dplyr::ungroup() %>% 
  dplyr::group_by(area_name, scenario_id, scenario_lab, cnstrnt_class) %>% 
  dplyr::mutate(
    patch_area_ha = as.numeric(sf::st_area(geometry))/10000
    , patch_area_rank = dplyr::dense_rank(dplyr::desc(patch_area_ha))
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::arrange(scenario_id, area_name, cnstrnt_class, dplyr::desc(patch_area_ha)) %>% 
  # calculate proportion of total area in each patch
  dplyr::inner_join(
    # total area of fireshed project areas included for analysis
      # after union to remove internal boundaries
    fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(scenario_id==1) %>% 
          dplyr::select(pa_id, area_name)
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      ) %>% 
      dplyr::group_by(area_name) %>% 
      dplyr::summarize(
        geometry = sf::st_union(geometry, by_feature = F)
      ) %>%
      dplyr::mutate(
        total_area_ha = as.numeric(sf::st_area(geometry))/10000
      ) %>% 
      sf::st_drop_geometry()
    , by = dplyr::join_by(area_name)
  ) %>% 
  dplyr::mutate(
    pct_patch_landscape = patch_area_ha / total_area_ha
  )
# function to highlight largest patch with fireshed project areas under
cols_rdylbu_dk <- c(
  "low constraint" = rev(RColorBrewer::brewer.pal(n=11,name="RdYlBu"))[2]
  , "med. constraint" = rev(RColorBrewer::brewer.pal(n=11,name="RdYlBu"))[7]
  , "high constraint" = rev(RColorBrewer::brewer.pal(n=11,name="RdYlBu"))[10]
)
largest_patch_map_fn <- function(nrow = 1){
    plt <- scenario_cnstrnt_map_fn(nrow) +
      geom_sf(
          data = cnstrnt_class_patches %>%
            dplyr::filter(
              area_name == area_nm_list[nrow]
              & patch_area_rank == 1
              & cnstrnt_class != "med. constraint"
            ) %>%
            sf::st_transform(crs=plt_crs)
          , aes(color = cnstrnt_class, fill = cnstrnt_class), lwd = 1.4
          , inherit.aes = F
          , alpha = 0
          # , show.legend = F
        ) +
      scale_color_manual(values = cols_rdylbu_dk, drop = F) +
      labs(
        color = "Level of\nConstraint"
        , subtitle = "largest patch of contiguous fireshed project areas classified by level of mechanical constraint" # largest patch of spatially adjacent...
      ) +
      theme(
        plot.subtitle = element_text(size = 8)
      ) +
      guides(
        color = "none"
        , fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9, linewidth = 0))
      )
  # return
  return(plt)
}
# largest_patch_map_fn(13)
# plot each landscape
plt_ptch_list <- 1:length(area_nm_list) %>%
  purrr::map(largest_patch_map_fn)
# scenario_maps
plt_ptch_list[]

Comparison of largest patch of lightly constrained fireshed project areas versus largest patch of highly constrained fireshed project areas.

Planning fewer, larger individual treatments across the landscape appears to be a better strategy when human community protection is a primary concern (Schmidt et al. 2008). Finney et al. (2007) note that the relative improvement of optimized treatment placement breaks down when larger proportions of the landscape (40–50%) are excluded from treatment because of land-management constraints.

# aggregate to constraint class level
cnstrnt_class_patches %>% 
  sf::st_drop_geometry() %>% 
  dplyr::group_by(area_name, scenario_lab, cnstrnt_class) %>% 
  dplyr::summarize(
    max_patch_area_ha = max(patch_area_ha, na.rm = T)
    , mean_patch_area_ha = mean(patch_area_ha, na.rm = T)
    , max_pct_patch_landscape = max(pct_patch_landscape, na.rm = T)
    , sum_pct_patch_landscape = sum(pct_patch_landscape)
  ) %>% 
  dplyr::ungroup() %>%
  dplyr::mutate(
    max_pct_patch_landscape_hack = dplyr::case_when(
      cnstrnt_class=="high constraint" ~ -max_pct_patch_landscape
      , TRUE ~ max_pct_patch_landscape
    )
  ) %>% 
  dplyr::filter(
    cnstrnt_class!="med. constraint"
  ) %>% 
# plot
ggplot(
  mapping = aes(
    y = scenario_lab
    , x = max_pct_patch_landscape_hack
    , fill = cnstrnt_class
    , group = cnstrnt_class
  )
) +
  # geom_line(
  #   mapping = aes(x = forcats::fct_rev(scenario_lab), y = max_pct_patch_landscape, color = forcats::fct_rev(cnstrnt_class), group = forcats::fct_rev(cnstrnt_class))
  #   , linewidth = 1.2
  # ) +
  geom_col(
    width = 0.7, alpha = 0.9
  ) +
  geom_vline(xintercept = 0, color = "gray55", linetype = "solid", linewidth = 1.1) +
  geom_text(
    mapping = aes(
      label = ifelse(
        max_pct_patch_landscape_hack<0
        , scales::percent(max_pct_patch_landscape, accuracy = 1)
        , ""
      )
      , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
      , fontface = "bold"
    )
    , color = "black", size = 2.7
    , hjust = +1
  ) +
  geom_text(
    mapping = aes(
      label = ifelse(
        max_pct_patch_landscape_hack>=0
        , scales::percent(max_pct_patch_landscape, accuracy = 1)
        , ""
      )
      , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
      , fontface = "bold"
    )
    , color = "black", size = 2.7
    , hjust = 0
  ) +
  annotate(
    geom = "text"
    , x = -0.77
    , y = 0
    , label = "high constraint"
    , color = cols_rdylbu_lt[3]
    , fontface = "bold"
    , size = 2.5
  ) +
  annotate(
    geom = "text"
    , x = 0.86
    , y = 0
    , label = "low constraint"
    , color = cols_rdylbu_lt[1]
    , fontface = "bold"
    , size = 2.5
  ) +
  scale_fill_manual(values = cols_rdylbu_lt) +
  scale_x_continuous(
    limits = c(-1.05,1.13)
    , labels = scales::percent_format()
    # , expand = expansion(mult = c(0, 0.1))
    # , breaks = scales::extended_breaks(5)
  ) +
  scale_y_discrete(expand = expansion(mult = c(0.8,0.3))) +
  facet_wrap(
    facets = vars(area_name)
    , ncol = 3
    , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    , scales = "free_y"
  ) +
  labs(
    fill = "Level of\nConstraint"
    , subtitle = "percentage of total landscape area comprised by the largest patch" # "Largest Patch % of Landscape Area"
    , x = "" # "Largest Patch % of Landscape Area"
    , y = ""
  ) +
  theme_light() +
  theme(
    legend.position = "none" # "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.ticks.x = element_blank()
    , panel.grid.major = element_blank()
    , panel.grid.minor = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 9)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

2.6 High-risk fireshed project areas

Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the plan (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).

2.6.1 Project Areas by WCS Status

# aggregate data by area, status
wcs_status_temp <- constrained_by_scnro_ls_pa %>% 
  dplyr::filter(scenario_id==1) %>%
  dplyr::group_by(area_name,fireshed_crisis_strategy) %>% 
  dplyr::summarise(
    n = n()
    , feature_area_ha = sum(feature_area_ha)
  ) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(
    tot_n = sum(n)
    , tot_area = sum(feature_area_ha)
    , pct_n = n/tot_n
    , pct_area = feature_area_ha/tot_area
  ) %>% 
  dplyr::ungroup()
# plot
ggplot(data = wcs_status_temp) +
  geom_col(
    mapping = aes(y = reorder(area_name,tot_n), x = n, fill = fireshed_crisis_strategy)
    , color = NA, width = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = area_name, x = n, group = fireshed_crisis_strategy
      , label = scales::comma(ifelse(n>20,n,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_pa %>% 
        dplyr::filter(scenario_id==1) %>%
        dplyr::group_by(area_name) %>% 
        dplyr::summarise(
          n = n()
          , feature_area_ha = sum(feature_area_ha)
        )
    , mapping = aes(
      y = reorder(area_name,n), x = n
      , label = scales::comma(n, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray")) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.1)),labels = scales::comma_format()) +
  labs(
    fill = "High-Risk\nStatus"
    , y = ""
    , x = "# of Fireshed\nProject Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
  )

2.6.2 Area of high-risk project areas

Only including fireshed project areas that have at least 25% of their area within the bounds of the priority landscapes and are in a high-risk fireshed.

See this section for discussion of high-risk firesheds and criteria for including fireshed project areas.

# data filter
wcs_status_temp %>% 
  dplyr::filter(fireshed_crisis_strategy!="Not High Risk") %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(tot_area = sum(feature_area_ha)) %>% 
# plot
ggplot() +
  geom_col(
    mapping = aes(y = reorder(area_name,tot_area), x = feature_area_ha, fill = fireshed_crisis_strategy)
    , color = NA, width = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = area_name, x = feature_area_ha, group = fireshed_crisis_strategy
      , label = ifelse(
          feature_area_ha<150000
          ,NA
          ,scales::comma(feature_area_ha,suffix = "k", scale = 1e-3, accuracy = 1)
        )
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_pa %>% 
        dplyr::filter(
          scenario_id==1
          & fireshed_crisis_strategy!="Not High Risk"
        ) %>%
        dplyr::group_by(area_name) %>% 
        dplyr::summarise(
          n = n()
          , feature_area_ha = sum(feature_area_ha)
        )
    , mapping = aes(
      y = reorder(area_name,feature_area_ha), x = feature_area_ha
      , label = scales::comma(feature_area_ha,suffix = "k ha", scale = 1e-3, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray")) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.15)),labels = scales::comma_format(suffix = "k ha", scale = 1e-3, accuracy = 1)) +
  labs(
    fill = "High-Risk\nStatus"
    , y = ""
    , x = "Area (ha) of Fireshed\nProject Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
  )

2.6.3 Reduction by Constraint

See this section for reduction by constraint at the landscape-level.

# aggregate by area_name, scenario for hi-risk pa's
constrained_by_scnro_ls_hionly_long <-
  constrained_by_scnro_ls_pa %>%
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::select(scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
                , (tidyselect::starts_with("rmn") & tidyselect::ends_with("area_ha"))
                , covertype_area_ha
  ) %>%
  dplyr::group_by(
    scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
  ) %>% 
  dplyr::summarise(
    dplyr::across(
      tidyselect::ends_with("area_ha")
      , sum
    )
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("rmn")
    , names_to = "constraint"
    , values_to = "rmn_area_ha"
    , names_prefix = "rmn"
    , values_drop_na = F
  ) %>%
  dplyr::mutate(constraint = stringr::str_remove(constraint,"_area_ha")) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , constraint = factor(
        constraint
        , ordered = TRUE
        , levels = c(
          "protected"
          , "slope"
          , "roads"
          , "riparian"
          , "administrative"
        )
        , labels = c(
          "Protected"
          , "Slope\nSteepness"
          , "Road\nDistance"
          , "Riparian\nBuffer"
          , "Administrative"
        )
      ) %>% forcats::fct_rev()
    , pct_rmn = rmn_area_ha/covertype_area_ha
    , pct_rdctn = dplyr::case_when(
      constraint_lvl == 1 ~ -1*(1-pct_rmn)
      , TRUE ~ -1*(dplyr::lag(pct_rmn)-pct_rmn)
    )
    , pct_rdctn_total = max(-1*(1 - ifelse(constraint_lvl == 5, pct_rmn,0)))
  )
  
# plot
ggplot(data = constrained_by_scnro_ls_hionly_long) +
  geom_col(
    mapping = aes(y = scenario_lab, x = pct_rdctn, fill = constraint)
    , color = NA, width = 0.7
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = scenario_lab, x = pct_rdctn
      , label = scales::percent(ifelse(pct_rdctn<.06*-1,pct_rdctn,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_hionly_long %>% dplyr::filter(constraint_lvl==5)
    , mapping = aes(
      y = scenario_lab, x = -1*(1 - pct_rmn)
      , label = scales::percent(-1*(1 - pct_rmn), accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_wrap(facets = vars(area_name)
      , ncol = 3
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    ) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0.02, .21)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
    , subtitle = "High-Risk Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 10) 
  )

The figure above shows the reduction by constraint for only high-risk fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape.

2.6.4 Reduction by Constraint high-risk vs. total

# combine data and plot
pct_rdctn_temp <- 
constrained_by_scnro_ls_hionly_long %>% 
  dplyr::ungroup() %>% 
  dplyr::select(scenario_lab, area_name, pct_rdctn, constraint_lvl) %>% 
  dplyr::mutate(dta_src = "High-Risk Project Areas") %>% 
  dplyr::bind_rows(
    constrained_by_scnro_ls_long %>% 
      dplyr::ungroup() %>% 
      dplyr::select(scenario_lab, area_name, pct_rdctn, constraint_lvl) %>% 
      dplyr::filter(!is.na(constraint_lvl)) %>% 
      dplyr::mutate(dta_src = "Overall Landscape Area")
  ) %>% 
  dplyr::group_by(scenario_lab, area_name, dta_src) %>% 
  dplyr::mutate(
    pct_rdctn_total = sum(pct_rdctn)
    , constraint = factor(
        constraint_lvl
        , ordered = TRUE
        , levels = 1:5
        , labels = c(
          "Protected"
          , "Slope\nSteepness"
          , "Road\nDistance"
          , "Riparian\nBuffer"
          , "Administrative"
        )
      ) %>% forcats::fct_rev()
    , scenario_lab = scenario_lab %>% forcats::fct_rev()
  )
# plot
ggplot(data = pct_rdctn_temp) +
  geom_col(
    mapping = aes(y = dta_src, x = pct_rdctn, fill = constraint)
    , color = NA, width = 0.7
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = dta_src, x = pct_rdctn
      , label = scales::percent(ifelse(pct_rdctn<.075*-1,pct_rdctn,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  geom_text(
    data = pct_rdctn_temp %>% dplyr::filter(constraint_lvl==5)
    , mapping = aes(
      y = dta_src, x = pct_rdctn_total
      , label = scales::percent(pct_rdctn_total, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_grid(
      cols = vars(scenario_lab)
      , rows = vars(area_name)
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
      , switch = "y"
    ) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0.02, .21)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_blank()
    , strip.text.x = element_text(color = "black", face = "bold")
    , strip.text.y.left = element_text(angle = 0)
    , strip.placement = "outside"
    
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

2.7 Objective vs. Treatable

Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the Strategy (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).

Using the area of fireshed project areas that have at least 25% of their area within the bounds of the priority landscapes and are in a high-risk fireshed as the treatment area objective (20-40%).

# Area of pa's by level of constraint
constrained_by_scnro_ls_cnstrnt_temp <- constrained_by_scnro_ls_pa %>% 
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::group_by(scenario_id,scenario_lab,state,name,area_name,cnstrnt_class) %>% 
  dplyr::summarise(
    n = n()
    , dplyr::across(
      c(feature_area_ha,covertype_area_ha,rmn5_administrative_area_ha)
      , ~ sum(.x,na.rm = T)
    )
  ) %>% 
  dplyr::group_by(scenario_id,scenario_lab,state,name,area_name) %>% 
  dplyr::mutate(
    # cnstrnt_class = stringr::word(cnstrnt_class) %>% stringr::str_remove_all("\\.")
    total_area_ha = sum(feature_area_ha)
    , total_covertype_area_ha = sum(covertype_area_ha)
    # "the need to treat 20 to 40 percent of the overall fireshed"
    , pct_treatable =  rmn5_administrative_area_ha/total_area_ha
    , objective_lo_treat_area_ha = total_area_ha*.2
    , objective_hi_treat_area_ha = total_area_ha*.4
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::arrange(state,name,scenario_id,cnstrnt_class)

2.7.1 Objective vs. Treatable Table

table_temp <-
  constrained_by_scnro_ls_cnstrnt_temp %>% 
  dplyr::mutate(
    cnstrnt_class = stringr::word(cnstrnt_class) %>% stringr::str_remove_all("\\.")
  ) %>% 
  dplyr::rename(
    remaining_pct=pct_treatable
    , remaining_ha=rmn5_administrative_area_ha
  ) %>%
  dplyr::select(-c(feature_area_ha, covertype_area_ha)) %>% 
  tidyr::pivot_wider(
    names_from = cnstrnt_class
    , values_from = c(n, remaining_pct, remaining_ha)
    , names_glue = "{cnstrnt_class}_{.value}"
    , values_fill = 0
  ) %>% 
  dplyr::mutate(
    tot_remaining_ha = low_remaining_ha+med_remaining_ha
    ,tot_remaining_pct = low_remaining_pct+med_remaining_pct
  ) %>% 
  dplyr::relocate(tot_remaining_ha,.after = med_remaining_ha) %>% 
  dplyr::inner_join(
    wf_landscapes %>% 
      sf::st_drop_geometry() %>% 
      dplyr::select(area_name,hectares) %>% 
      dplyr::rename(landscape_area_ha=hectares)
    , by = dplyr::join_by(area_name)
  )
# make table 
kableExtra::kable(
    table_temp %>% 
      dplyr::select(
        c(
          scenario_lab
          # , landscape_area_ha
          , total_area_ha
          , total_covertype_area_ha
          # light and mod area
          , low_remaining_ha, med_remaining_ha, tot_remaining_ha
          , low_remaining_pct, med_remaining_pct, tot_remaining_pct
        )
      ) %>% 
      dplyr::mutate(
        dplyr::across(
          tidyselect::ends_with("_ha")
          , ~ scales::comma(.x, suffix = "k", scale = 1e-3, accuracy = 1)
        )
        , dplyr::across(
          tidyselect::ends_with("_pct")
          , ~ scales::percent(.x, accuracy = .1)
        )
      )
  , caption = "Wildfire Crisis Strategy High-Risk Fireshed Project Areas<br>mechanically available for treatment in lightly and moderately constrained project areas (≥25% in landscape)"
  , col.names = c(
    ""
    # , "Total\nLandscape"
    , "Total Area (ha)"
    , "Forest+Shrub Area (ha)"
    , "Lightly Constrained", "Moderately Constrained"
    , "Total"
    , "Lightly Constrained", "Moderately Constrained"
    , "Total"
  
  )
  , escape = F
) %>%  
    add_header_above(c(" " = 1, "High-Risk Fireshed\nProject Areas" = 2, "Treatable High-Risk\nArea (ha)"=3, "Treatable High-Risk\nPercent (%)" = 3)) %>%
    kable_classic(full_width=T) %>% 
    pack_rows(index = table(forcats::fct_inorder(table_temp$area_name))) %>% 
    kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%  
    kableExtra::column_spec(
      9 # 10
      , color = "white"
      , background = 
          table_temp %>% 
            dplyr::mutate(
              is_gt20 = dplyr::case_when(
                low_remaining_pct >= .2 ~ "#008b8b"
                , tot_remaining_pct >= .2 ~ "#545454"
                , TRUE ~ "#cd3700"
              )
          ) %>% 
          dplyr::pull(is_gt20)
          
    ) %>% 
    kableExtra::scroll_box(width = "740px")
Table 2.3: Wildfire Crisis Strategy High-Risk Fireshed Project Areas
mechanically available for treatment in lightly and moderately constrained project areas (≥25% in landscape)
High-Risk Fireshed
Project Areas
Treatable High-Risk
Area (ha)
Treatable High-Risk
Percent (%)
Total Area (ha) Forest+Shrub Area (ha) Lightly Constrained Moderately Constrained Total Lightly Constrained Moderately Constrained Total
AZ: 4FRI
Scenario 1 1,579k 1,436k 360k 85k 445k 22.8% 5.4% 28.2%
Scenario 2 1,579k 1,436k 554k 58k 611k 35.0% 3.7% 38.7%
Scenario 3 1,579k 1,436k 936k 30k 966k 59.3% 1.9% 61.2%
AZ: Prescott
Scenario 1 221k 191k 34k 19k 53k 15.3% 8.6% 23.9%
Scenario 2 221k 191k 69k 14k 83k 31.3% 6.2% 37.6%
Scenario 3 221k 191k 112k 0k 112k 50.9% 0.0% 50.9%
AZ: San Carlos Apache Tribal Forest Protection
Scenario 1 250k 230k 4k 29k 33k 1.6% 11.6% 13.2%
Scenario 2 250k 230k 42k 36k 78k 16.8% 14.4% 31.2%
Scenario 3 250k 230k 80k 20k 100k 31.9% 7.9% 39.8%
CA: North Yuba
Scenario 1 206k 197k 32k 26k 59k 15.8% 12.8% 28.6%
Scenario 2 206k 197k 85k 6k 91k 41.4% 2.8% 44.2%
Scenario 3 206k 197k 104k 7k 111k 50.5% 3.6% 54.1%
CA: Plumas Community Protection
Scenario 1 157k 148k 24k 25k 49k 15.3% 15.7% 31.0%
Scenario 2 157k 148k 68k 5k 72k 43.1% 3.0% 46.1%
Scenario 3 157k 148k 91k 0k 91k 57.9% 0.0% 57.9%
CA: Southern California Fireshed Risk Reduction Strategy
Scenario 1 1,053k 854k 44k 80k 124k 4.2% 7.6% 11.8%
Scenario 2 1,053k 854k 186k 68k 255k 17.7% 6.5% 24.2%
Scenario 3 1,053k 854k 253k 65k 318k 24.0% 6.1% 30.2%
CA: Stanislaus
Scenario 1 153k 143k 31k 21k 52k 20.4% 13.7% 34.1%
Scenario 2 153k 143k 63k 6k 70k 41.4% 4.1% 45.5%
Scenario 3 153k 143k 84k 0k 84k 54.9% 0.0% 54.9%
CA: Trinity Forest Health and Fire-Resilient Rural Communities
Scenario 1 394k 336k 0k 30k 30k 0.0% 7.5% 7.5%
Scenario 2 394k 336k 32k 37k 68k 8.1% 9.3% 17.4%
Scenario 3 394k 336k 140k 20k 160k 35.6% 5.0% 40.6%
CO: Colorado Front Range
Scenario 1 1,138k 823k 120k 100k 220k 10.5% 8.8% 19.3%
Scenario 2 1,138k 823k 301k 57k 358k 26.5% 5.0% 31.4%
Scenario 3 1,138k 823k 398k 43k 441k 34.9% 3.8% 38.7%
ID: Nez Perce-Clearwater-Lower Salmon
Scenario 1 520k 444k 91k 39k 130k 17.4% 7.6% 25.0%
Scenario 2 520k 444k 203k 12k 215k 39.1% 2.2% 41.3%
Scenario 3 520k 444k 213k 8k 221k 40.9% 1.6% 42.5%
ID: Southwest Idaho
Scenario 1 734k 632k 243k 32k 275k 33.1% 4.3% 37.4%
Scenario 2 734k 632k 408k 20k 428k 55.6% 2.7% 58.3%
Scenario 3 734k 632k 429k 11k 441k 58.4% 1.5% 60.0%
MT: Kootenai Complex
Scenario 1 322k 294k 101k 13k 114k 31.5% 4.0% 35.5%
Scenario 2 322k 294k 143k 6k 149k 44.5% 1.9% 46.4%
Scenario 3 322k 294k 194k 7k 200k 60.2% 2.0% 62.3%
NV: Sierra and Elko Fronts
Scenario 1 678k 453k 119k 39k 157k 17.5% 5.7% 23.2%
Scenario 2 678k 453k 231k 17k 248k 34.1% 2.5% 36.7%
Scenario 3 678k 453k 257k 13k 270k 37.9% 2.0% 39.9%
NM: Enchanted Circle
Scenario 1 318k 285k 51k 29k 79k 16.0% 9.0% 25.0%
Scenario 2 318k 285k 122k 13k 135k 38.3% 4.1% 42.5%
Scenario 3 318k 285k 142k 6k 149k 44.7% 2.0% 46.7%
OR: Central Oregon
Scenario 1 278k 232k 145k 5k 150k 52.2% 1.8% 54.0%
Scenario 2 278k 232k 169k 2k 170k 60.5% 0.6% 61.2%
Scenario 3 278k 232k 203k 0k 203k 73.1% 0.0% 73.1%
OR: Klamath River Basin
Scenario 1 971k 808k 33k 63k 96k 3.4% 6.5% 9.8%
Scenario 2 971k 808k 190k 47k 237k 19.5% 4.9% 24.4%
Scenario 3 971k 808k 316k 48k 364k 32.5% 4.9% 37.4%
OR: Mount Hood Forest Health and Fire-Resilient Communities
Scenario 1 339k 244k 22k 22k 44k 6.4% 6.6% 13.0%
Scenario 2 339k 244k 49k 21k 71k 14.6% 6.3% 20.9%
Scenario 3 339k 244k 135k 8k 143k 39.8% 2.3% 42.1%
UT: Pine Valley
Scenario 1 187k 173k 18k 20k 39k 9.9% 10.9% 20.8%
Scenario 2 187k 173k 53k 12k 65k 28.5% 6.3% 34.9%
Scenario 3 187k 173k 57k 13k 69k 30.3% 6.8% 37.2%
UT: Wasatch
Scenario 1 156k 126k 23k 10k 33k 14.8% 6.6% 21.4%
Scenario 2 156k 126k 45k 8k 52k 28.6% 5.0% 33.7%
Scenario 3 156k 126k 50k 6k 56k 31.9% 4.0% 36.0%
WA: Central Washington Initiative
Scenario 1 818k 588k 21k 32k 52k 2.5% 3.9% 6.4%
Scenario 2 818k 588k 76k 44k 120k 9.3% 5.3% 14.6%
Scenario 3 818k 588k 258k 29k 287k 31.6% 3.6% 35.1%
WA: Colville Northeast Washington Vision
Scenario 1 272k 244k 70k 13k 83k 25.9% 4.7% 30.6%
Scenario 2 272k 244k 122k 8k 129k 44.7% 2.9% 47.6%
Scenario 3 272k 244k 149k 6k 155k 54.7% 2.2% 56.9%

2.7.2 Reduction in treatable area flow

table_temp %>% 
  dplyr::filter(scenario_id==1) %>% 
  dplyr::select(
    c(
      area_name
      , landscape_area_ha
      , total_area_ha
      , total_covertype_area_ha
      # , tot_remaining_ha
      # , low_remaining_ha
    )
  ) %>%
  tidyr::pivot_longer(
    !area_name
    , names_to = "name"
    , values_to = "hectares"
  ) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(order = dplyr::row_number(),y_val=2) %>% 
  dplyr::ungroup() %>% 
  dplyr::cross_join(data.frame(scenario_id = 1:3)) %>% 
  ####
  dplyr::bind_rows(
    table_temp %>% 
      dplyr::select(
        c(
          area_name
          , scenario_id
          , tot_remaining_ha
          , low_remaining_ha
        )
      ) %>%
      tidyr::pivot_longer(
        !c(area_name,scenario_id)
        , names_to = "name"
        , values_to = "hectares"
      ) %>% 
      dplyr::group_by(area_name,scenario_id) %>% 
      dplyr::mutate(
        order = dplyr::row_number()+3
        , y_val = scenario_id
      ) %>% 
      dplyr::ungroup() %>% 
      dplyr::relocate(scenario_id,.after = dplyr::last_col())
  ) %>% 
  dplyr::arrange(area_name,order,y_val ) %>% 
  dplyr::mutate(
    order_lab = factor(
      order
      , levels = 1:5
      , labels = c(
        "Total\nLandscape"
        , "High-risk\nFireshed PAs"
        , "High-risk\nForest+Shrub"
        , "Lt.+Md. Constrained"
        , "Lightly Constrained"
      )
      , ordered = T
    )
  ) %>% 
  dplyr::inner_join(
    wf_landscapes %>% 
      sf::st_drop_geometry() %>% 
      dplyr::select(area_name, hectares) %>% 
      dplyr::rename(tot_hectares=hectares)
    , by = dplyr::join_by(area_name)
  ) %>% 
  dplyr::inner_join(
    table_temp %>% 
      dplyr::filter(scenario_id==1) %>% 
      dplyr::select(area_name, objective_lo_treat_area_ha)
    , by = dplyr::join_by(area_name)
  ) %>% 
  dplyr::mutate(
    y_val = -y_val
    , sizing = ifelse(hectares>tot_hectares,1,hectares/tot_hectares)
    , coloring = dplyr::case_when(
        order %in% 4:5
          & hectares >= objective_lo_treat_area_ha
          ~ 1
        , order %in% 4:5
          & hectares < objective_lo_treat_area_ha
          ~ 2
        , T ~ 3
    )
  ) %>% 
  # dplyr::filter(stringr::str_starts(area_name,"CA")) %>% 
  # View()
ggplot() +
  geom_line(
    mapping = aes(x = order_lab, y = y_val, group = scenario_id)
    , color = "gray"
  ) +
  geom_point(
    mapping = aes(x = order_lab, y = y_val, size = sizing*100, color = factor(coloring))
    , shape = 15
  ) +
  geom_text(
    mapping = aes(
      x = order_lab, y = y_val, group = scenario_id
      , label = scales::comma(hectares,suffix = "k", scale = 1e-3, accuracy = 1)
    )
    , size = 2.5
    , color = "black"
    , vjust = 2.5
    , angle = 90
  ) +
  geom_text(
    mapping = aes(
      x = order_lab, y = y_val, group = scenario_id
      , label = ifelse(
        order == 4, paste0("Scenario", y_val*-1), NA
      )
      , fontface = "italic"
    )
    , size = 2.5
    , color = "black"
    , hjust = 0.1
    , vjust = -1.4
  ) +
  facet_wrap(facets = vars(area_name)
    , ncol = 3
    , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    , scales = "free_x"
  ) +
  scale_color_manual(values = c("#008b8b", "#cd3700", "#bebebe")) +
  scale_y_continuous(expand = expansion(c(0.3,0.3))) +
  scale_x_discrete(labels = scales::label_wrap(10)) +
  labs(
    x = ""
    , y = ""
    , subtitle = "<span><span style='color:#008b8b;'><b><i>≥20% treatable</i></b></span> | <span style='color:#cd3700;'><b><i><20% treatable</i></b></span></span>"
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.y = element_blank()
    , axis.ticks.y = element_blank()
    , axis.text = element_text(size = 7)
    , panel.grid = element_blank()
    , plot.subtitle = ggtext::element_markdown(size = 10)
    , strip.text = element_text(color = "black", face = "bold", size = 10)
  )

2.7.3 Objective vs. Treatable Plot

obj_trtbl_temp <-
  constrained_by_scnro_ls_cnstrnt_temp %>% 
  dplyr::filter(cnstrnt_class!="high constraint") %>% 
  dplyr::left_join(
    # aggregate by area_name, scenario_lab
    constrained_by_scnro_ls_cnstrnt_temp %>% 
        dplyr::filter(cnstrnt_class!="high constraint") %>% 
        dplyr::group_by(area_name, scenario_id) %>% 
        dplyr::mutate(low=ifelse(cnstrnt_class=="low constraint",pct_treatable,0)) %>% 
        dplyr::summarise(
          pct_treatable_lowmed = sum(pct_treatable)
          , pct_treatable_low = sum(low)
        ) %>% 
        dplyr::mutate(
          treatment_objective = dplyr::case_when(
              pct_treatable_low >= .2 ~ 1
              , pct_treatable_lowmed >= .2 ~ 2
              , TRUE ~ 3
            ) %>% 
            factor(
              levels = 1:3
              , labels = c("mechanical alone","mechanical+fire","fire")
              , ordered = T
            )
        ) %>% 
        dplyr::group_by(scenario_id) %>% 
        dplyr::arrange(scenario_id, treatment_objective
                       , desc(pct_treatable_low), desc(pct_treatable_lowmed)
        ) %>% 
        dplyr::mutate(sorting_rank = dplyr::row_number()) %>% 
        dplyr::ungroup()
      , by = dplyr::join_by(area_name, scenario_id)
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(
    dplyr::across(
      c(cnstrnt_class, scenario_lab)
      , ~ forcats::fct_rev(.x)
    )
  ) 
# plot
obj_plt_fn <- function(scnum){
plt <- ggplot(
    data = obj_trtbl_temp %>% dplyr::filter(scenario_id==scnum)
    , mapping = aes(
        x = pct_treatable
        , y = reorder(area_name, desc(sorting_rank))
      )
  ) +
  geom_rect(mapping = aes(xmin = 0.2, xmax=0.4, ymin = -Inf, ymax = Inf)
    , fill = "gray77"
    , alpha = 0.6
  ) +
  annotate("text"
     , x = 0.46
     , y = obj_trtbl_temp %>% 
          dplyr::filter(scenario_id == scnum & sorting_rank==max(obj_trtbl_temp$sorting_rank)) %>% 
          dplyr::pull(area_name) %>% 
          unique()
     , label = expression(bold("objective")~bold("20-40%"))
     , size = 3, color = "gray65", parse = T
    ) +
  geom_col(
    mapping = aes(fill = cnstrnt_class)
    , color = NA, width = 0.7
  ) +
  geom_text(
    mapping = aes(
      group = cnstrnt_class
      , label = scales::percent(
          ifelse(pct_treatable<.06,NA,pct_treatable)
          , accuracy = 1
        )
      # , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.2
  ) +
  # total
  geom_text(
    data = obj_trtbl_temp %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::group_by(area_name, scenario_id) %>% 
        dplyr::filter(dplyr::row_number()==1) %>% 
        dplyr::ungroup()
    , mapping = aes(
      x = pct_treatable_lowmed
      , y = reorder(area_name, desc(sorting_rank))
      , label = scales::percent(pct_treatable_lowmed, accuracy = 1) 
      , color = treatment_objective
      , fontface = "bold"
    )
    , size = 3.5
    , hjust = -0.1
    , show.legend = F
  ) +
  # facet_grid(
  #   rows = vars(scenario_lab)
  #   , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
  #   , scales = "free_y"
  #   # , switch = "y"
  # ) +
  scale_fill_manual(values = RColorBrewer::brewer.pal(n=3,name="RdYlBu")[2:3]) +
  scale_color_manual(values = c(
    "mechanical alone"="darkcyan"
    ,"mechanical+fire"="gray33"
    , "fire"="orangered3")
  ) +
  # scale_x_continuous(expand = expansion(mult = c(0, .15)),labels = scales::percent_format()) +
  scale_x_continuous(limits = c(0, .78),labels = scales::percent_format(),position = "top") +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "Percent of Total High-Risk Area"
    , subtitle = obj_trtbl_temp %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::pull(scenario_lab) %>% unique()
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=8)
    , axis.text.x = element_blank()
    , axis.ticks.x = element_blank()
    , axis.text = element_text(size = 8)
    , panel.grid.major.x = element_blank()
    , panel.grid.minor.x = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 10)
    , strip.background = element_rect(fill = "gray88")
    # , strip.placement = "outside"
    , plot.margin = margin(0, 0, 0, 0, "cm")
    , plot.subtitle = element_text(vjust = -4,face="bold")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )
return(plt)
}
cowplot::plot_grid(
  obj_plt_fn(1) +
    theme(
      legend.position = c(0.86,0.82)
      , legend.direction = "vertical"
    )
  , obj_plt_fn(2) +
    theme(legend.position = "none")
  , obj_plt_fn(3) +
    # scale_x_continuous(position = "top") +
    # labs(x = "Percent of Total High-Risk Area") +
    theme(legend.position = "none", plot.margin = margin(0, 0, 0.1, 0, "cm"))
  , ncol = 1
  , align = 'vh'
)

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_obj_sc.png")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

2.8 Totals

##################
# Data Prep
##################
# union landscape and pa level data
# aggregate by scenario for hi-risk pa's
pct_rdctn_pas_temp <- 
  constrained_by_scnro_ls_pa %>%
  dplyr::group_by(pa_id,scenario_id) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::select(scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
                , (tidyselect::starts_with("rmn") & tidyselect::ends_with("area_ha"))
                , covertype_area_ha, feature_area_ha
  ) %>%
  dplyr::group_by(scenario_id, scenario_desc, scenario_lab) %>% 
  dplyr::summarise(
    dplyr::across(
      tidyselect::ends_with("area_ha")
      , sum
    )
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("rmn")
    , names_to = "constraint"
    , values_to = "rmn_area_ha"
    , names_prefix = "rmn"
    , values_drop_na = F
  ) %>%
  dplyr::mutate(constraint = stringr::str_remove(constraint,"_area_ha")) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , pct_rmn = rmn_area_ha/covertype_area_ha
    , pct_rdctn = dplyr::case_when(
      constraint_lvl == 1 ~ -1*(1-pct_rmn)
      , TRUE ~ -1*(dplyr::lag(pct_rmn)-pct_rmn)
    )
    , pct_rdctn_total = max(-1*(1 - ifelse(constraint_lvl == 5, pct_rmn,0)))
  ) %>% 
  dplyr::ungroup()
# overall
pct_rdctn_ls_temp <- 
  constrained_by_scnro_ls %>%
  dplyr::select(scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
                , (tidyselect::starts_with("rmn") & tidyselect::ends_with("area_ha"))
                , covertype_area_ha, feature_area_ha
  ) %>%
  dplyr::group_by(scenario_id, scenario_desc, scenario_lab) %>% 
  dplyr::summarise(
    dplyr::across(
      tidyselect::ends_with("area_ha")
      , sum
    )
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("rmn")
    , names_to = "constraint"
    , values_to = "rmn_area_ha"
    , names_prefix = "rmn"
    , values_drop_na = F
  ) %>%
  dplyr::mutate(constraint = stringr::str_remove(constraint,"_area_ha")) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , pct_rmn = rmn_area_ha/covertype_area_ha
    , pct_rdctn = dplyr::case_when(
      constraint_lvl == 1 ~ -1*(1-pct_rmn)
      , TRUE ~ -1*(dplyr::lag(pct_rmn)-pct_rmn)
    )
    , pct_rdctn_total = max(-1*(1 - ifelse(constraint_lvl == 5, pct_rmn,0)))
  ) %>% 
  dplyr::ungroup()
# combine
pct_rdctn_all_temp <- 
  pct_rdctn_pas_temp %>% 
  dplyr::mutate(dta_src = "High-Risk Project Areas") %>% 
    dplyr::bind_rows(
      pct_rdctn_ls_temp %>% 
        dplyr::mutate(dta_src = "Overall Landscape Area")
    ) %>% 
    dplyr::mutate(
      constraint_lab = factor(
          constraint_lvl
          , ordered = TRUE
          , 1:5
          , labels = c(
            "Protected"
            , "Slope\nSteepness"
            , "Road\nDistance"
            , "Riparian\nBuffer"
            , "Administrative"
          )
        ) %>% forcats::fct_rev()
    ) %>% 
    dplyr::relocate(dta_src)

2.8.1 Reduction Treatable Area Table

tot_table_temp <- pct_rdctn_all_temp %>% 
  dplyr::mutate(constraint_nm = paste0(constraint_lvl,constraint)) %>% 
  dplyr::select(-c(constraint_lvl,constraint,constraint_lab)) %>% 
  tidyr::pivot_wider(
    names_from = constraint_nm
    , values_from = c(rmn_area_ha, pct_rmn, pct_rdctn)
  ) %>% 
  dplyr::mutate(
    pct_covertype_area = covertype_area_ha/feature_area_ha
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::contains("area_ha")
      , ~ scales::comma(.x,suffix = " M", scale = 1e-6, accuracy = 0.1)
    )
    , dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ scales::percent(.x, accuracy = 0.1)
    )
  ) %>% 
  dplyr::arrange(
    desc(dta_src)
    , scenario_id
  )
# make table 
  kableExtra::kable(
      tot_table_temp %>% 
        dplyr::select(
          scenario_lab
          , covertype_area_ha
          , pct_covertype_area
          , pct_rdctn_1protected
          , pct_rdctn_2slope
          , pct_rdctn_3roads
          , pct_rdctn_4riparian
          , pct_rdctn_5administrative
          , rmn_area_ha_5administrative
          , pct_rmn_5administrative
        )
      , caption = paste0("<b>Wildfire Crisis Strategy Priority Landscapes and High-Risk Fireshed Project Areas</b><br>percent reduction of different types of constraints on mechanical treatment")
      , col.names = c(
        ""
        , "Forest+Shrub\n(ha)", "Forest+Shrub\n%"
        , "Protected"
        , "Slope\nSteepness"
        , "Road\nDistance"
        , "Riparian\nBuffer"
        , "Administrative"
        , "Remaining (ha)"
        , "Remaining (%)"
      )
      , escape = F
    ) %>%  
    add_header_above(c(" " = 1, "Strategy Priority\nArea"=2, "Constraints on\nMechanical Treatment" = 5, " " = 2)) %>%
    kable_classic(full_width=T) %>% 
    pack_rows(index = table(forcats::fct_inorder(tot_table_temp$dta_src))) %>% 
    kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%  
    kableExtra::scroll_box(width = "740px")
Table 2.4: Wildfire Crisis Strategy Priority Landscapes and High-Risk Fireshed Project Areas
percent reduction of different types of constraints on mechanical treatment
Strategy Priority
Area
Constraints on
Mechanical Treatment
Forest+Shrub (ha) Forest+Shrub % Protected Slope Steepness Road Distance Riparian Buffer Administrative Remaining (ha) Remaining (%)
Overall Landscape Area
Scenario 1 15.8 M 83.1% -21.9% -16.9% -19.9% -7.0% -5.5% 4.6 M 28.9%
Scenario 2 15.8 M 83.1% -21.9% -5.2% -12.1% -10.5% -8.2% 6.7 M 42.2%
Scenario 3 15.8 M 83.1% -21.9% -5.2% -12.1% -7.5% 0.0% 8.4 M 53.4%
High-Risk Project Areas
Scenario 1 8.7 M 82.5% -17.5% -21.1% -18.2% -8.5% -6.3% 2.5 M 28.5%
Scenario 2 8.7 M 82.5% -17.5% -6.4% -10.5% -12.9% -9.4% 3.7 M 43.2%
Scenario 3 8.7 M 82.5% -17.5% -6.4% -10.5% -9.3% 0.0% 4.9 M 56.3%

2.8.2 Reduction by Constraint

# plot
# plot
ggplot(data = pct_rdctn_all_temp) +
  geom_col(
    mapping = aes(y = scenario_lab, x = pct_rdctn, fill = constraint_lab)
    , color = NA, width = 0.6
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = scenario_lab, x = pct_rdctn
      , label = scales::percent(ifelse(pct_rdctn<.06*-1,pct_rdctn,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = pct_rdctn_all_temp %>% dplyr::filter(constraint_lvl==5)
    , mapping = aes(
      y = scenario_lab, x = pct_rdctn_total
      , label = scales::percent(pct_rdctn_total, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_grid(
      cols = vars(dta_src %>% forcats::fct_rev())
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    ) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0, .12)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_blank()
    , strip.text.x = element_text(color = "black", face = "bold")
    , strip.placement = "outside"
    
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

2.8.3 Distribution of Fireshed PA Constraint

cnstrnt_dta_temp <- constrained_by_scnro_ls_pa %>% 
  dplyr::group_by(pa_id,scenario_id) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::filter(fireshed_crisis_strategy  %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::count(scenario_id, scenario_desc, scenario_lab, cnstrnt_class) %>% 
  dplyr::group_by(scenario_desc, scenario_lab) %>% 
  dplyr::mutate(
    pct = n/sum(n)
    , dta_src = "High-Risk Firesheds"
  ) %>% 
  dplyr::bind_rows(
    constrained_by_scnro_ls_pa %>% 
      dplyr::group_by(pa_id,scenario_id) %>% 
      dplyr::filter(dplyr::row_number()==1) %>% 
      dplyr::ungroup() %>% 
      dplyr::count(scenario_id, scenario_desc, scenario_lab, cnstrnt_class) %>% 
      dplyr::group_by(scenario_desc, scenario_lab) %>% 
      dplyr::mutate(
        pct = n/sum(n)
        , dta_src = "Overall Landscape Area"
      )
  ) %>% 
  dplyr::mutate(
    dta_src = dta_src %>% forcats::fct_rev()
  )
# n for notation
nlab_temp <- 
  cnstrnt_dta_temp %>% 
  dplyr::filter(scenario_id==1) %>% 
  dplyr::group_by(dta_src) %>% 
  dplyr::summarise(n = sum(n)) %>% 
  dplyr::mutate(l = paste0(dta_src, " n = ", scales::comma(n))) %>% 
  dplyr::pull(l) %>% 
  paste(collapse = "\n")
# plot base
plt_fshed_cnstrnt_lvl_temp <- ggplot(data=cnstrnt_dta_temp) +
  geom_col(
    mapping = aes(x = pct, y = scenario_lab, fill=cnstrnt_class)
    ,width = 0.7, alpha=0.8
  ) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  facet_grid(cols = vars(dta_src)) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "% of Fireshed\nProject Areas"
    # , caption = nlab_temp
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
    , plot.caption = element_text(size = 8)
  ) + 
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )
# plot with large labels
plt_fshed_cnstrnt_lvl_temp + 
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.065,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 4 # 2.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = ifelse(pct>=0.065,
          paste0("n="
            ,scales::comma(n, accuracy = 1)
          )
          , NA)
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.5 # 2.8
    , vjust = 2.3
  ) 

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_prjareas_dist_sc.png")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

2.8.3.1 Combine map and fireshed plot

# plot small for inset
inset_temp <- plt_fshed_cnstrnt_lvl_temp +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.105,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3.2 # 2.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = ifelse(pct>=0.105,
          paste0("n="
            ,scales::comma(n, accuracy = 1)
          )
          , NA)
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.1 # 2.8
    , vjust = 2.3
  ) +
    cowplot::theme_minimal_grid(9) +
    theme(
      legend.position = "none"
      , plot.background = element_rect(color = "black", fill="white", linewidth = 2)
      , axis.text.x = element_blank()
      , axis.title.x = element_text(size=7)
      , strip.text = element_text(color = "black", face = "bold", size = 9)
    )

# combine
cowplot::ggdraw(
  plt_fpas_sc_map[[2]] + 
    theme(
      legend.position = c(.4,.98)
      , legend.direction = "horizontal"
      , legend.text = element_text(size = 9)
      , plot.background = element_rect(fill="white", color = NA)
      , plot.margin = margin(0, 5.1, -0.7, 0.1, "cm")
    )
  ) +
  cowplot::draw_plot(inset_temp, x=.5, y=.56, width=.45, height=.4) +
  cowplot::draw_plot_label(
    label = c("A", "B"),
    x = c(0.02, 0.5),
    y = c(1, 0.95),
    size = 17,
    fontface = "bold",
    color = "gray11"
  )

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_prjareas_sc_map_inset.png")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

2.9 Landcover Class Area

2.9.1 Overall Landscape

What is the landcover classification of the area within the landscape boundary? This GEE script was used to calculate the area by NLCD landcover class.

# nlcd class groups
nlcd_reclass_df_temp <- data.frame(
    class_id = c(11,12,21,22,23,24,31,41,42,43,51,52,71,72,73,74,81,82,90,95) %>% as.character()
    , class_label = c(
      "Water"
      , "Water"
      , "Developed"
      , "Developed"
      , "Developed"
      , "Developed"
      , "Barren"
      , "Forest"
      , "Forest"
      , "Forest"
      , "Shrubland"
      , "Shrubland"
      , "Herbaceous"
      , "Herbaceous"
      , "Herbaceous"
      , "Herbaceous"
      , "Planted/Cultivated" # double-check name
      , "Planted/Cultivated"
      , "Wetlands"
      , "Wetlands"
      )
  ) %>% 
  dplyr::left_join(
    FedData::nlcd_colors() %>% dplyr::rename(class_id=ID) %>% dplyr::rename_with(tolower) %>% dplyr::mutate_all(as.character)
    , by = dplyr::join_by(class_id)
  )
# load data by wf priority 
wf_nlcd_area <- readr::read_csv("../data/wildfirepriority/wfpriority_nlcd_area.csv") %>% 
  dplyr::rename_with(tolower) %>% 
  dplyr::left_join(
    wf_landscapes %>% 
      dplyr::select(name,state,area_name)
    , by = dplyr::join_by(name,state)
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("area_m2_nlcd_cl")
    , names_to = "class_id"
    , values_to = "area_m2"
    , names_prefix = "area_m2_nlcd_cl_"
    , values_drop_na = F
  ) %>% 
  dplyr::left_join(
    nlcd_reclass_df_temp
    , by = dplyr::join_by(class_id)
  ) %>% 
  dplyr::mutate(
    class_label = forcats::fct_relevel(class_label, c("Forest","Shrubland")) %>% forcats::fct_rev()
  )
# get color palette
nlcd_colors <- wf_nlcd_area %>% 
  dplyr::group_by(class_id,color,class_label) %>% 
  dplyr::summarise(sum_area = sum(area_m2,na.rm = T)) %>% 
  dplyr::arrange(class_label,desc(sum_area)) %>% 
  dplyr::group_by(class_label) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::pull(color)
# aggregate and plot
wf_nlcd_area %>% 
  dplyr::group_by(area_name,class_label) %>% 
  dplyr::summarise(sum_area_ha = sum(area_m2,na.rm = F)/10000) %>% 
  dplyr::mutate(sum_area_ha = ifelse(is.na(sum_area_ha),0,sum_area_ha)) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(
    pct_area = sum_area_ha/sum(sum_area_ha)
    # , area_name = area_name %>% forcats::fct_rev()
  ) %>%
  dplyr::ungroup() %>% 
ggplot(
    mapping = aes(x = pct_area, y = reorder(area_name, desc(area_name)), fill=class_label)
  ) +
  geom_col(width = 0.7, alpha=0.8) +
  geom_text(
    mapping = aes(
      label = scales::percent(ifelse(pct_area>.03,pct_area,NA), accuracy = 1)
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  scale_fill_manual(values = nlcd_colors) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "NLCD\nCover\nType"
    , y = ""
    , x = "% Landscape Area"
  ) +
  theme_light() +
  theme(
    legend.position = "top" # c(0.9, 0.9)
    , legend.title = element_text(size=9)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

2.9.2 Fireshed Project Areas

What is the landcover classification of the fireshed project area units? This GEE script was used to calculate the area by NLCD landcover class.

# load nlcd cover data by fireshed project area
nlcd_area_ls_pa <- 
  readr::read_csv("../data/wildfirepriority/wfpriority_fireshed_nlcd_area.csv") %>% 
  dplyr::rename_with(tolower) %>% 
  dplyr::select(
    pa_id
    , tidyselect::contains("area_m2")
  ) %>% 
  dplyr::mutate(pa_id = as.character(pa_id)) %>% 
  dplyr::inner_join(
    constrained_by_scnro_ls_pa %>% 
      dplyr::filter(scenario_id==1) %>% 
      dplyr::select(pa_id,name,state,area_name,fireshed_crisis_strategy)
    , by = dplyr::join_by(pa_id)
    , multiple = "all"
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("area_m2_nlcd_cl")
    , names_to = "class_id"
    , values_to = "area_m2"
    , names_prefix = "area_m2_nlcd_cl_"
    , values_drop_na = F
  ) %>% 
  dplyr::left_join(
    nlcd_reclass_df_temp
    , by = dplyr::join_by(class_id)
  ) %>% 
  dplyr::mutate(
    class_label = forcats::fct_relevel(class_label, c("Forest","Shrubland")) %>% forcats::fct_rev()
  )
# aggregate and plot
nlcd_area_ls_pa %>% 
  dplyr::group_by(area_name,class_label) %>% 
  dplyr::summarise(sum_area_ha = sum(area_m2,na.rm = F)/10000) %>% 
  dplyr::mutate(sum_area_ha = ifelse(is.na(sum_area_ha),0,sum_area_ha)) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(pct_area = sum_area_ha/sum(sum_area_ha)) %>%
  dplyr::ungroup() %>% 
ggplot(
    mapping = aes(x = pct_area, y = reorder(area_name, desc(area_name)), fill=class_label)
  ) +
  geom_col(width = 0.7, alpha=0.8) +
  geom_text(
    mapping = aes(
      label = scales::percent(ifelse(pct_area>.03,pct_area,NA), accuracy = 1)
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  scale_fill_manual(values = nlcd_colors) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "NLCD\nCover\nType"
    , y = ""
    , x = "% Landscape Area"
    , subtitle = "All Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top" # c(0.9, 0.9)
    , legend.title = element_text(size=9)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

The figure above shows the NLCD landcover class distribution for all fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape.

2.9.3 NLCD Cover Type high-risk firesheds

# aggregate and plot
nlcd_area_ls_pa %>%
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::group_by(area_name,class_label) %>% 
  dplyr::summarise(sum_area_ha = sum(area_m2,na.rm = F)/10000) %>% 
  dplyr::mutate(sum_area_ha = ifelse(is.na(sum_area_ha),0,sum_area_ha)) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(pct_area = sum_area_ha/sum(sum_area_ha)) %>%
  dplyr::ungroup() %>% 
ggplot(
    mapping = aes(x = pct_area, y = reorder(area_name, desc(area_name)), fill=class_label)
  ) +
  geom_col(width = 0.7, alpha=0.8) +
  geom_text(
    mapping = aes(
      label = scales::percent(ifelse(pct_area>.03,pct_area,NA), accuracy = 1)
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  scale_fill_manual(values = nlcd_colors) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "NLCD\nCover\nType"
    , y = ""
    , x = "% Landscape Area"
    , subtitle = "High-Risk Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top" # c(0.9, 0.9)
    , legend.title = element_text(size=9)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

The figure above shows the NLCD landcover class distribution for only high-risk fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape.

2.9.4 NLCD Cover Type vs. Mechanically Available

Is there a relationship between forest and shrubland cover and mechanical availability and is that relationship the same?

nlcd_area_ls_pa_cl_temp <-
  nlcd_area_ls_pa %>% 
    # dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
    dplyr::group_by(pa_id,area_name,class_label) %>% 
    dplyr::summarise(sum_area_ha = sum(area_m2,na.rm = F)/10000) %>% 
    dplyr::group_by(pa_id,area_name) %>% 
    dplyr::mutate(pct_area = sum_area_ha/sum(sum_area_ha)) %>%
    dplyr::ungroup() %>% 
    dplyr::filter(class_label %in% c("Forest","Shrubland")) %>% 
    dplyr::inner_join(
      constrained_by_scnro_ls_pa %>% 
        dplyr::select(
          scenario_id
          , scenario_lab
          , scenario_desc
          , pa_id
          , area_name
          , pct_rmn5_administrative
          , pct_rdctn_total
          , cnstrnt_class
          , fireshed_crisis_strategy
        )
      , by = dplyr::join_by(pa_id,area_name)
      , multiple = "all"
    )  
  # dplyr::mutate(scenario_desc = scenario_desc %>% forcats::fct_rev())
# treatable vs covertype
ggplot(
    data = nlcd_area_ls_pa_cl_temp
    , mapping = aes(y = pct_rmn5_administrative, x = pct_area, color = class_label)
  ) +
  geom_point(alpha = 0.7, size = 1) +
  # geom_smooth(method = "lm", lwd = 1, color = "gray", linetype = "dashed", se=F) +
  scale_color_manual(values = nlcd_colors[(length(nlcd_colors)-1):length(nlcd_colors)]) +
  facet_grid(cols = vars(class_label), rows = vars(scenario_desc)) +
  scale_x_continuous(labels = scales::percent_format()) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    x = "Percent Cover Type"
    , y = "Percent Treatable Area Remaining"
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , axis.title = element_text(size=9)
    , axis.text = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
  )

All fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape shown in figure above.

2.9.5 Distribution of Mechanically Available (Shrb vs. Frst)

# n for notation
nlab_temp <- nlcd_area_ls_pa_cl_temp %>% 
  dplyr::filter(pct_area > .5 & scenario_id==1) %>% 
  dplyr::count(class_label) %>% 
  dplyr::mutate(l = paste0(class_label, " n = ", scales::comma(n))) %>% 
  dplyr::pull(l) %>% 
  paste(collapse = " | ")
# plot
nlcd_area_ls_pa_cl_temp %>% 
  dplyr::filter(pct_area > .5) %>% 
ggplot(
    mapping = aes(x = pct_rmn5_administrative, fill = class_label, group = class_label, color = class_label)
  ) +
  geom_density(alpha = 0.7) +
  scale_fill_manual(values = nlcd_colors[(length(nlcd_colors)-1):length(nlcd_colors)]) +
  scale_color_manual(values = nlcd_colors[(length(nlcd_colors)-1):length(nlcd_colors)]) +
  facet_grid(cols = vars(class_label), rows = vars(scenario_desc)) +
  scale_x_continuous(labels = scales::percent_format()) +
  # scale_y_continuous(labels = scales::percent_format()) +
  labs(
    y = "Density"
    , x = "Percent Treatable Area Remaining"
    , caption = nlab_temp
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , axis.title = element_text(size=9)
    , axis.text = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
    , plot.caption = element_text(size = 8)
  )

In the plot above fireshed project areas are grouped into forest or shrubland based on the majority cover type (e.g. >50% forest). Fireshed project areas where the majority covertype is not forest or shrubland are not shown. All fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape shown in figure above.

2.9.6 Distribution of Fireshed PA Constraint (Shrb vs. Frst)

nlcd_area_ls_pa_cl_temp %>% 
  dplyr::filter(pct_area > .5) %>% 
  dplyr::count(class_label, scenario_desc, scenario_lab, cnstrnt_class) %>% 
  dplyr::group_by(class_label, scenario_desc, scenario_lab) %>% 
  dplyr::mutate(pct = n/sum(n)) %>% 
ggplot() +
  geom_col(
    mapping = aes(x = pct, y = scenario_lab, fill=cnstrnt_class)
    ,width = 0.7, alpha=0.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.05,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  facet_grid(cols = vars(class_label)) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "% of Fireshed\nProject Areas"
    , caption = nlab_temp
    , subtitle = "All Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
    , plot.caption = element_text(size = 8)
  ) + 
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

In the plot above fireshed project areas are grouped into forest or shrubland based on the majority cover type (e.g. >50% forest). Fireshed project areas where the majority covertype is not forest or shrubland are not shown. All fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape shown in figure above.

2.9.7 High-Risk Distribution of Fireshed PA Constraint (Shrb vs. Frst)

# n for notation
nlab_temp <- nlcd_area_ls_pa_cl_temp %>% 
  dplyr::filter(pct_area > .5 
                & fireshed_crisis_strategy  %in% c("USFS-Only", "All-Lands") 
                & scenario_id==1) %>% 
  dplyr::count(class_label) %>% 
  dplyr::mutate(l = paste0(class_label, " n = ", scales::comma(n))) %>% 
  dplyr::pull(l) %>% 
  paste(collapse = " | ")
# plot
nlcd_area_ls_pa_cl_temp %>% 
  dplyr::filter(pct_area > .5 & fireshed_crisis_strategy  %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::count(class_label, scenario_desc, scenario_lab, cnstrnt_class) %>% 
  dplyr::group_by(class_label, scenario_desc, scenario_lab) %>% 
  dplyr::mutate(pct = n/sum(n)) %>% 
ggplot() +
  geom_col(
    mapping = aes(x = pct, y = scenario_lab, fill=cnstrnt_class)
    ,width = 0.7, alpha=0.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.05,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  facet_grid(cols = vars(class_label)) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "% of Fireshed\nProject Areas"
    , caption = nlab_temp
    , subtitle = "High-Risk Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
    , plot.caption = element_text(size = 8)
  ) + 
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

In the plot above fireshed project areas are grouped into forest or shrubland based on the majority cover type (e.g. >50% forest). Fireshed project areas where the majority covertype is not forest or shrubland are not shown. Only high-risk fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape shown in figure above.